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ABSTRACT 


In order to solve FD-TD electromagnetic radiation and scattering problems in free 
space, the computational grid should be truncated at some finite distance to simulate the 
solution for an infinitely large grid. One way to emulate an infinite grid is to absorb 
outgoing waves incident onto the grid boundary such that there is no reflection back into 
the grid. Waves incident to the boundary must "exit" the grid boundary without affecting 
the solution within the grid. Based on this idea, this thesis develops Transparent Grid 
Termination (TGT) in 1-D, 2-D and 3-D as an absorbing boundary. TGT performance is 
compared with a very large grid, for which reflections do not return to the computational 
domain within the observation time window. 


v 





TABLE OF CONTENTS 


l. INTRODUCTION . 1 

A. OVERVIEW . 1 

B. OBJECTIVE . 3 

H. ANALYSIS OF TGT FOR 1-D FD-TD . 7 

A. FD-TD FORMULATION IN 1-D . 7 

1. FD-TD Equations in 1-D . 7 

2. Derivation of Magnetic Field Update Equation . 10 

3. Derivation of Electric Field Update Equation . 15 

B. TRANSPARENT GRID TERMINATION IN 1-D . 20 

1. 1-D Grid . 20 

2. Grid Termination . 22 

C. 1-D TGT RESULTS . 28 

m. ANALYSIS OF TGT FOR 2-D FDTD . 31 

A. FD-TD FORMULATION IN 2-D . 31 

1. FD-TD Equations in 2-D . 31 

2. Derivation of TEz Magnetic Field Update Equation . 35 

3. Derivation of TEz Electric Field Update Equations . 40 

4. Derivation of TMz Electric Field Update Equation . 47 

5. Derivation of TMz Magnetic Field Update Equations . 53 

B. TRANSPARENT GRID TERMINATION IN 2-D . 58 

1. 2-D TEz and TMz Grids . 58 

2. 2-D Grid Termination . 62 

C. 2-D TGT RESULTS . 74 

IV. ANALYSIS OF TGT FOR 3-D FD-TD . 91 

A FD-TD FORMULATION IN 3-D .. 91 

1. FD-TD Equations in 3-D . 91 

2. Derivation of Magnetic Field Update Equations . 94 

3. Derivation of Electric Field Update Equations . 100 































B. TRANSPARENT GRID TERMINATION IN 3-D . 108 

1. 3-D Grid . 108 

2. 3-D Grid Termination. 114 

C. 3-D TGT RESULTS . 124 

V. SUMMARY AND CONCLUSION . 129 

A. SUMMARY . 129 

B. CONCLUSIONS AND RECOMMENDATIONS . 129 

APPENDIX . 131 

A. DBIR.M . 131 

B. TGT TM.M . 132 

LIST OF REFERENCES . 135 

INITIAL DISTRIBUTION LIST . 137 


viii 















I. INTRODUCTION 


A. OVERVIEW 

The analysis of electromagnetic problems in the time domain has become more 
common in recent decades. In this approach, the differential or integral equations for the 
particular problem under consideration are typically solved numerically using a short 
time-duration waveform as excitation. By far the most common of these techniques is the 
so-called finite difference time-domain (FD-TD) method. 

In the FD-TD method, the problem is spatially discretized ("sampled") using a spatial 
step Al, usually in a Cartesian coordinate system and temporally "sampled" using a time 
step At. The differential operators in the differential form of Maxwell's equations are 
approximated by finit e differences or, equivalently, integrals in the integral form of 
Maxwell's equations are approximated by finite sums. In this research, a pulse waveform is 
excited at the center of the discretization grid, and the fields at the grid nodes are 
computed at discrete time steps. The field at a particular grid node at time t can be 
determined from the fields of this and adjacent nodes calculated for the previous time step 
t - At. The electric and magnetic fields obtained via the FD-TD method are the result of a 
"marching in time" process whose accuracy is affected by the approximations introduced 
at each time step. These include the grid size, time step, length of time observation, 
excitation waveform, and so forth. 
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There are several advantages to the FD-TD method. First, it is relatively simple to 
implement for complex objects because of spatial discretization. Inhomogeneous media is 
handled by assigning varying electric and magnetic properties (s, p, and a) to individual 
cells. Although the FD-TD electric and magnetic field "update" equations appear complex, 
they contain only addition, substraction, and multiplication. Furthermore, computer 
memory requirements are usually less demanding than those for the method of moments 
(MOM). 

An important consideration in the FD-TD method is how to terminate the spatial 
discretization grid. The FD-TD equation derived in the following section will describe 
how waves propagate along an infinit e grid. However, in all practical applications the grid 
will have to be restricted to a finite number of nodes. Therefore, a method for grid 
"termination" is required. The "standard" FD-TD update equations are valid for all nodes 
except those on the grid edges, because the grid edge nodes have fewer neighbors than the 
non-edge (internal) nodes. Therefore, edge nodes with a different number of neighbors 
have update equations different from the equations derived for the internal grid nodes. The 
question is how can one dete rmin e the equations for the edge nodes? One approach is to 
try to absorb a wave incident from inside the grid onto the boundary such that there is no 
reflection back into the grid. This is referred to as an absorbing boundary condition 
(ABC) [Ref 2], However, the grid itself is not ideal, that is a wave does not pass through 
the grid without some reflection even for those nodes that are not at grid's edges (this is 
the consequence of spatial and temporal field sampling). Requiring that there is no 
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reflection off the grid edges is thus too restrictive, because there is some reflection at each 
and every node within the grid. We thus propose a different approach that we refer to as 
the Transparent Grid Termination (TGT). We start by making only one demand on the 
grid termination condition: that, ideally, the solution within the grid should be the same 
solution that one would obtain with an infinitely large grid, or in other words, that 
reflections off a grid edge node are the same as the reflections off any other node inside 
the grid. Note that if the grid were infinitely large, the grid edges would be at infinity and 
the equations for the grid edge nodes would be irrelevant, since the excited or scattered 
field would take an infinitely long time to reach the grid edges, and thus the edge 
equations would never really be needed. The TGT requirement is conceptually very simple 
and fairly straightforward to implement. We will use the multiport concept (used in 
analysis of linear systems) to introduce the TGT concept, and details will be shown in the 
following sections. 

B. OBJECTIVE 

The performance of the TGT in 1-D, 2-D, and 3-D can be assessed based on the 
residual reflection off the TGT boundary. In this thesis the total power of the 
electromagnetic field, excited by a unit amplitude delta pulse in the grid center, will be 
observed as a function of time. If the grid were ideal, the energy of the field at any node 
will be zero for all time after the pulse has passed that node. The total energy within an 
infinit e ideal grid would be constant however. The total energy within a finite ideal grid 
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with ideal grid termination condition will be constant until the pulse has reached the grid 
edge and will then quickly go to zero as the pulse leaves the grid. However, because the 
grid is not ideal, there are some residual values at any grid node, even after the initial pulse 
has passed that node. Therefore, the energy within a non-ideal grid will not be constant 
but will fluctuate about its value for the ideal grid. These fluctuations generally decrease 
with time but they would exist even if the grid were infinite. If the grid is terminated, the 
residual energy within the grid will increase because of additional reflections off the grid 
edges. The global effect of the grid termination can thus be assessed by comparing the 
total grid powers as functions of time for the following two cases: an "infinite" grid, and a 
finite size grid with TGT or ABC. Note that the "infinite" grid is simulated by using a 
large grid and finishin g the observations before the pulse that started at the grid's cotter 
has reached any of the grid edge nodes. Note that if the TGT or ABC were ideal, there 
would be no increase in energy within the grid after the pulse has reached the grid edges. 
If the TGT is not ideal, the relative change of the energy within the grid (compared to the 
case of an " infinit e" grid) is a measure of TGT "quality": the smaller the change, the better 
the TGT would be. We will thus compare the total energies within the TGT boundary for 
a large grid with Dirichlet boundary condition (any boundary condition is appropriate for 
this " infinit e 11 grid simulation because the reflections off the edges of the large grid do not 
reach the smaller grid boundary within the observation time-window) and the smaller grid 
with the TGT boundary. Once again, the large grid must provide sufficient distance 
between the TGT and Dirichlet boundaries such that the reflections of the Dirichlet 
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boundary do not reach the TGT boundary within the observation interval. The measure of 
residual reflection off the TGT boundary will be the difference of the energies within the 
TGT boundary for the " infinit e 11 grid and the TGT. It is convenient to normalize this 
difference to the energy within the TGT boundary for the "infinite" grid, and express the 
resulting ratio in dB. The objective of this thesis is thus to determine the quality 
(measured in dB as described above) of the Transparent Grid Termination (TGT). 
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H. ANALYSIS OF TGT FOR 1-D FD-TD 


A. FD-TD FORMULATION IN 1-D 

1. FD-TD Equations in 1-D 

The incident and the scattered electromagnetic fields and the media parameters in 

1-D problems can vary with one spatial coordinate only. We will denote this coordinate 

as z. The fields will also be functions of time t. The media will be assumed stationary, 

that is, the media parameters are not functions of time t. 1-D electromagnetic fields must 

— > — > 

be Transverse Electro-Magnetic (TEM) fields, that is the field components £ and if 
must be in a plane transverse to the direction of field propagation [Ref 5], Since we have 

denoted the direction of propagation as z, the fields can be denoted as transverse-to-z or 

TEM Z . The electric and magnetic vectors of a TEM field must be perpendicular to each 

other. The E and H field vectors and the direction of propagation form a triplet of 

orthogonal vectors, as shown below. 
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Figure 2.1. TEM Field Components. 

The axes were selected such that the electric field vector is x-directed and the 
magnetic field vector is y-directed. 

E(z,t)=E,(z,f)~x (2.1) 

H(z,f)=H,(z,t)y (2.2) 

The electric and magnetic fields satisfy Maxwell's curl equations which express the 

electric and magnetic field coupling and can be written in the integral or the differential 

form. We will use the integral forms because they are applicable to finite (line, surface, or 

volume) domains whereas the differential forms are applicable to infinitesimally small 

domains (points). Shown below are the integral forms of Maxwell's curl equations for a 

TEM Z field. 

E x (z, t)~x • dl E = -J~ J]^ \x(z)H y (z, t)~y • ds E ■ (2.3) 

| c H y (z, i)~y • dl H = ^jjj^ s (z)E x (z, t)~x • ds H ■ +JJ^ c s(z)E x (z, t jx • ds H (2.4) 
The second equation does not have the source current term since it has beea 
assumed that there were no source currents in the domain of interest. The contours of 
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integration for the electric and the magnetic field contour integrals are, in general, different 
and are thus labeled C E for an electric field circulation and C H for a magnetic field 
circulation. (A closed path line integral will be referred to as the circulation.) Similarly, the 
surfaces associated with the contours are labeled S E for the magnetic field surface integral 
and S H for the electric field surface integral. Surface integrals on the right-hand side will be 
referred to as fluxes. The first step in discretizing the curl equations is to select the 
contours for the electric and magnetic field circulation. Since the essence of Maxwell's 
curl equations is the coupling of the electric and magnetic fields, the contours will be 
selected such that this coupling is achieved in a straightforward manner, as shown below. 



Figure 2.2. Contours for Electric and Magnetic Field Circulation and Fluxes. 

The electric field contours C E and the magnetic field contours C H are square (A1 by 
Al) contours, but in orthogonal planes. The C E contours are in the xz-plane, while the C H 
contours are in the yz-plane. The contours can be compared to link s of a chain, evoking 
the idea of electric and magnetic field linkage. This selection of such contours leads to the 
so-called dual discretization grid, also referred to as the Yee lattice [Ref. 1]. 
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2. Derivation of Magnetic Field Update Equation 


The discrete equivalent of the electric field circulation will be determined next. A 
contour C E is shown below. 


x 

A 


/ “ — 

E(z-Al/2,t) * H(z,t) E(z-+Al/2,t) 


Figure 2.3. A Contour for Electric Field Circulation and Magnetic Flux Calculations. 

The center of the surface S E is assumed to have the coordinates (x, z). A local 
coordinate system (£,Q can be established, with the origin at the center of the contour. 
Any point (x',z?) within or on the contour C E (or, equivalently, any point on the surface 

S E ) can be specified by its local coordinates B, and C, 

x‘ -x + B, z' = z + C, where and 

The local coordinates will be used in evaluation of the line and surface integrals 

that constitute the integral foims of Maxwell's equations. The electric field E x is not a 

function of the x-coordinate and the electric field is zero along the top and the bottom 

sides of the contour because there can be no z-directed field components of a TEM Z field. 

The circulation of the electric field around the electric field contour C E , assuming 

counter-clockwise reference direction such that the normal to the surface S E is in the 

—y 

direction of the magnetic field H(z, t), can be evaluated exactly 
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| E x (z +C, t)x *dt E = [*l E x (z + f, t)x • d£x -\ + l E x (z- f, t)x • (2.5) 

*2 2 

E x (z+^t)x .~dt E = [E x (z+f,t)-E x (z-f,t)]-Al (2.6) 

The rate of change of magnetic flux through the surface S E (within the contour C E ) 
can, on the other hand, be evaluated only approximately, since the exact way that the 
magnetic flux density 

B(z, t) = \i(z)H(z, t) (2.7) 

varies with the z-coordinate is not known a priori. Note that there is no variation of the 

magnetic flux density with the x-coordinate. We will next evaluate the magnetic flux 

approximately 

IT n(z+C)^+C,()7-3 I = f! (2-8) 

JJSg 7 7 

The surface integral reduces to a line integral, since the integrand is not a function 
of the local x-coordinate S, 

\i(z + QH y (z + C, t)~y *2 E =Al £* p(z+Qff,(z+C,')<*C (2.9) 

There are infinit ely many ways to model the variation of the magnetic flux density 
with the local z-coordinate C, within the contour. Since the magnetic flux density B(z+£,t) 
is the product of the permeability p(z+Q and the magnetic field H(z+^,t) we first need to 
assume a certain variation of the magnetic field with C, such that the integral can be 
evaluated over S E . The simplest model assumes that the contour width A1 is small enough 
such that the magnetic field H(z+£,t) within the contour may be assumed constant and 
equal to the value at the contour's center H(z,t). This is equivalent to using a piece-wise 
constant ("pulse" expansion) approximation of the actual magnetic field variation with the 
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z-coordinate. The above assumption allows that the magnetic field, although constant 
within a contour, can change from one contour to an other. This yields an approximate 
expression for the for the magnetic flux 

(f (2.10) 

Sg 2 

The integral of the permeability p(z+Q can be re-written in the following manner 

f| n( Z + M = A/-[i|| rtz+C«] (2.11) 

The term in the brackets is recognized as the average permeability (i^ within the contour 

C E - 

H^(r) = iJ^|n(2+g< (212) 

The approximate expression for the magnetic flux through S E can now be written using the 
average permeability as 

ff (i(z + C)Hy {~ + C, f) y • dsp; ~ (A/) \iavg(z)H y (z,i) (2-13) 

Sg 

Again, this approximate expression resulted from the piece-wise constant 
approximation of the magnetic field with respect to the z-coordinate. The main advantage 
of the piece-wise constant expansion employed above is its simplicity. Better accuracy can 
be achieved by using more involved models for the field variation with z, but at the 
expense of increasing the complexity and computational time. The first curl equation of 

Maxwell can now be replaced by an approximate equation 

[E,(z+f,t)-E,(z- f.0] • M 

which can be simplified (because of media stationarity) to 

E x (z + y, f) ~ E x (z — y, /) « {H y (z, ?)} • A/ • llavgi^) (2.15) 

This equation can also be re-wntten as 
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( 2 . 16 ) 


The time derivative operator in the above equation is typically replaced by the 
finite difference approximation [Ref 1]. However, an alternate approach is presented here, 
such that the approximation of the field temporal variation is shown to be analogous to the 
approximations already introduced for the field spatial variation. The above equation is 
integrated with respect to time to get an approximate equation for the magnetic field at the 
present time t 

Hy(z ’ 0 * A/- pl g (a) [Jo Ex( ~ z + a ’ ■ ' 11 Ex(z ~ * < 2 ' 17 > 

A similar integral equation can be written for the magnetic field at the same spatial 
location but at an earlier time t-At 

A' ■- AO « + 7 • - r EJt-- f. X)A] (2.18) 

Subtracting the two equations we get 

H y (z,t)-H y (z,t-At)* E x (z + f, t )<h-E x (z- f, x)ch] (2.19) 

which can be also written as an "update" equation (assuming that the previous field value 
at the same location is known) 

H„(z,t) « Hy(z,I- 41)- a/ .^ (z) [L&<?+ f• -0*-L-f.Tjrfr] (2.20) 
The integrals on the right-hand side can not be evaluated exactly, because the 
exact temporal variation of the electric fields within the At interval prior to t is generally 
not known, just like the spatial integral for the magnetic flux could not have been 
evaluated exactly because the exact variation of the magnetic field within the contour was 
not known. However, the integrals can be evaluated approximately by assuming a certain 
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variation of the electric field with the temporal variable t. The simplest approach, 
consistent with the assumptions made for the field spatial variation, would be to assume 
that the interval At is small enough such that the electric field can be assumed 
approximately constant within At and equal to the value at the center of the time interval 
(t-At, t) 


E(z + f,z)■■ 

:E<Z + f,t-f) 

for 

t-At <x <t 

(2.21) 

E(Z- f.T)« 


for 

t-At<x<t 

(2.22) 


The above represents a piece-wise constant approximation of the electric field variation 
with respect to the temporal variable t. The approximate expression for the magnetic field 


at location z and the time t now becomes 


H,(z, t) « Hy(z, t-At)- + f, t - f) - E,<? - f, t - f)] (2.23) 

The above equation can be written in a somewhat different form, using the identity 

1 ... 1 1 ^ i 1 1 1_L = v „Ii = v xi 0 24) 

V VoVr ys °Z 0 Vr H' ' 


=-p=7T = v o4~lr = v oTotJ- (2.24) 

,n8n /Lin M-r Zo V* M-r 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Zo « 3770), and p r denotes relative permitivity. Introducing the grid 
"propagation" velocity 


_ A / 

V grid ~ Af 


(2.25) 


we can write the approximate "update" equation for the magnetic field as 

H r (z, t) * Hy(z, t-At)- +f . t-- f) ■- E,(z- f, t-- f)] (2.26) 

The equation simplifies for the case of non-magnetic media (p=l) 

H y (z, t) * H y (z, t-At)- Yo^[E x (z + f,t-f)-E x (z-f,t-f)] (2.27) 
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The last two equations show that the magnetic field H and the electric field E are 
evaluated at points shifted spatially by Al/2, and at instants separated in time by At/2. The 
relationship between spatial and temporal "samples" of the electric and the magnetic fields 
is thus the same the samples are shifted with respect to each other by one-half of the 
sampling interval. This is the essence of the Yee lattice and applies to "sampling" of 2-D 
and 3-D fields as well. Dual approximate equations for the electric field "updates" can 
now be derived using the same procedure shown for the magnetic field update equations. 

3. Derivation of Electric Field Update Equation 

We start with the Maxwell's curl equation for the circulation of the magnetic field 
j Cif H y (z, t)~x • dl H = ~ s (z)E x (z, t)~x • ds H + JJ^ o(z)E x (z, t)~? • ds H (2.28) 

The discrete equivalent of the magnetic field circulation will be determined first, 
using the contour shown below. 

Ch 

1 / .—.> z 

I * H(z-tAl/2,t) 


Figure 2.4. A Contour for Magnetic Field Circulation and Electric Flux Calculations. 

A local coordinate system (v|/,Q can be established, with the origin at the center of 
the contour. Any point (y\z') within or on the contour C H (or, equivalently, any point on 

the surface S H ) can be specified by its local coordinates 

y' =y + vj/ z / =z + 4 where -^<v|/<+^ and - ^ < C - 







The circulation of the magnetic field around C H , assuming counter-clockwise 

reference direction such that the normal to the surface S H is in the direction of the electric 
—^ 

field E (z, t), is therefore given exactly by the following simple expression 

I H y (z + &)7'dh=\i,Hy(z-7,<)7'‘lvy (2-29) 

JC H *2 2 

I H y (z + &tf •dh = [H y (z-f,t)-Hy(z + f,t)]-Al (2.30) 

* '-'H 

The rate of change of the electric flux through the surface S H can only be evaluated 

approximately, since the exact way that the electric flux density 

D(z,t) = z(z)H(z,t ) (2.31) 

varies in the z-direction is not known a priori. First we evaluate the flux integrals 

ff s(z + QExiz + t)~x t f + I s(z + QE x (z + C, t)~x . ~xd\\)d^ (2.32) 

JJ Stf 2 2 

ff <j(z + QE x (z +C, t)~x • dsn = \ + l o(z+QE x (z + Q, t)~x • ~xd\\/d £ (2.33) 

JJ Sff 2 ~2 

Since the integrands are not functions of the local y-coordinate vj/ the surface integrals 
simplify to line integrals 

ff e(z + QE x (z+ C, t)~x • dsn = A/ f l s(z + QE x (z + Q, t)dC, (2.34) 

JJSfi J ~2 

ff a (z + QE x (z + C, t)~x •ds„ = Al f + £ o(z + QE x (z + £, t)rf£ (2.35) 

Next, we assume that the contour width A1 is small enough such that the electric 
field E(z+£,t) within the contour may be assumed constant and equal to the value at the 
contour's center E(z,t). This assumption yields approximate expressions for the for the 
flux integrals 

ff s(z + QE x (z + C, t)~x • ds H = A/ • E x (z, t) • fj} s (z+QdC (2.36) 
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(2.37) 


jj SK c(z+ QE x (z +C t)x • ~ds „= M■ E,(z, I) • a(z+CK 

or, using the average permitivity and conductivity 

jj^ s(z+ QE x (z +C ,• ds H = {Atf - Zav g (z) • £*(z, 0 (2.38) 

j}^ ar(z + QE x (z + 0^ •<&// = (A/) 2 • o<»s( z ) * E x {z, t) (2.39) 

The second curl equation of Maxwell can now be replaced by an approximate equation 
\H y {z - f, t)-H y {z + f, 0] • A/ * |{ [s^z^z, 0](A0 2 } + (A/) 2 • a^(z) • £,(z, t) (2.40) 
which can be simplified (because of media stationarity) to 

H y {z - f, t)-H y (z + f,/) « £<n>g(z) • A/ • |{£,(z, t)} + A/• c^(z) • E x (z, t) (2.41) 


The approximate equation can be re-written as 

- t. *>■-*# + * •')] - si £ ' (z ’ 0 <242) 

The above equation is integrated with respect to time to get an approximate equation for 
the electric field at the present time t 

E - (! ’ 0 * H ' {! -T-W- f. "'< z + 7. H f, E < 243 > 

A .similar equation can be written for the electric field at the same spatial location 
but at an earlier time t-At 

< 2 - 44 > 

Subtracting the two equations we get 

E x {z,t)*E x (z,t-At)+ (2.45) 

1 


Al-Savgiz) 


;l - f. ■*>* -- L H ><? + f • ■»>*] -- Si L £ ' <z> T)rfI 
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Assuming that the interval At is small enough such that the magnetic fields can be 
assumed approximately constant within At and equal to the value at the center of the time 
interval At 


we get 


H(z + f,x)> 

*H(z + f,t-f) for t-At<x<t 

(2.46) 

H(z- 

*H(z-f,t-f ) for t - At <x <t 

(2.47) 


E x (z,t)xE x (z,t-At)+ 

(2.48) 


A t r u /_ A/ , Af\ U {. I A/ , Af\l Gqygi 2 ) f' 

+ ,)J Sm(z) i, 

The above equation can be written in a somewhat different form 

E x (z, t) « E x (z, t - At)+ 

At r tu /_ A/ . A/\ rr , A/ ( Ar\"l A /Tavg( Z )r 1 




(2.49) 


The term in the brackets is recognized as the average value of the electric field 
within the interval (t-At, t). The above equation can be written in a more compact form 


using the identity 


1 _ 1 1 _ 7^ 1 1 _ 1 /Ho I z I 

e eos r JixtJeEJZZSr V £ o 0 0£ ' 


(2.50) 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Zo « 3770) and e t denotes relative permitivity. Introducing the grid 
"propagation" velocity 

* % (2 - 51) 

we can write the approximate equation for the electric field "updates" as 

E x (z,t)*E x (z,t-At)+ (2.52) 


vp 1 

V grid g r aV g(z) 


H y (z -f,t- f) - H y {z + j,t-j)-Al-Gavg(z) j- f L E x (z, x)dX- 
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or 

E x (z,t)*E x (z,t-At)+ (2.53) 

Zo^~ ~Hy( z + T’*~ y) - A/ - a avg (z)E Xy avg(^ t )(z, /)J 

The average of the electric field between t-At and t, denoted E x>lv ^ t) can not be 
calculated exactly because the exact temporal variation of the electric field is not known a 
priori. However, assuming a linear variation within At interval, the average field is the 


average of the values 


ai rat 




Wife vv idhlzAZlEtM 


(2.54) 


; E_,. t) (z) into the update equation we get 


Substituting the time-average E x _ avg(At) ( 

E x (z, t) « E x (z, t - A t)+ (2.55) 

<7 V 0 1 EJ /_ A/ # A/\ rj , 4/ # A t\ * j s \ E x {z,t — Af) E x {z,f)~\ 

Zo ^ - 7.y) - ^+T>' - t) - -5-J 

which gives the final equation for the electric field updates as 

1 1 vp Zo • A/ • Oqvg(z) 

< 2 - 56 ) 

J + 1 v 0 

2 ^grid £ r,avg (z) 

2 q v o 1 

v grid Srfivg(z) fur. A/ , Af\ u /_ , A/ * Af\”| 

1 vo 2o-A/-q w # g ' (, -7’ < -T)- g ^- | -T.<-7)J 
+ 2 £r,av^(z) 

The equation simplifies greatly for the case of non-conductive media (a=0) 

E,(z, 0 * E,(z, I - Ai) + Zo^- j7 -^[h,(z -¥,/-§)- H y (z + f, r- f)] (2.57) 

In the case of free-space (s=l), the equation simplifies further 

Z x (z, t - At) + Zo ^[H y (z - f, / •- f )- H y {z + f, / - f)] (2.58) 


E x (zJ) «E X 
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These electric field "update" equations also show that the electric field E and the 
magnetic field H should be evaluated at points shifted spatially by Al/2, and at instants 
separated in time by At/2. 

B. TRANSPARENT GRID TERMINATION IN 1-D 

1. I-D Grid 

The electric and magnetic field update equations have been derived using local 
coordinate systems, with origins at the centers of the magnetic and electric contours, 
respectively. These equations now need to be "converted" to a global coordinate system, 
that is to the grid of equi-distant sampling points along the z-axis (for 1-D). We will 
assume that our domain ("grid") is a line segment of length L. The fields are sampled 
using a spatial step A1=L/N z . The electric and magnetic field sample locations are 
"interleaved", as shown below for N z =5. The first and the last spatial sampling points 
form the grid edges. The spatial edge samples can be either electric field samples or 
magnetic field samples. Although the selection of the field for the grid edges makes no 
difference in principle, we will in general use electric field samples for the grid edges. 

#.O.•.O.•.O.•.O.•.o.• 

A1 

i«.> 

Figure 2.5. 1-D E-field (black) and H-field (white) nodes. 
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We would evaluate the fields within the grid first, since a known incident field 


propagates from the center. Electric field update equations can now be "converted" to 

electric field grid equations by replacing the variables z and t with the grid coordinates of 

the electric field spatial and temporal sampling points 

z -» JcAl, k = 0,1 ...N z and t nAt , n = 0,1...N, 

The electric field update equation for non-conductive media is 

E x (kAl, nAi) « E x (kAl, nAt - At) (2.59) 

* ■ nA " » ■ 

The notation can be simplified further by omitting the common A1 and At terms, and using 


a superscript for the index of the temporal sampling point 


£;(*)»£r'(*)+ZoQ£-) 


1 


v grids Z navg (k) 


H n ;\k-\)-H7~ l {k + \) 


(2.60) 


The electric field grid equation for conductive media, using the same notation as above, is 


E”{k) 


1 if Vo ^ Z 0 * A/ ’ <5avg(k) 

1 “ 2yvgridJ' 


•r,a\g 


(k) 


1 , if v 0 ^ Zo-Al-Oavgik) 

2 yVgridJ S r ,avg(k) 


ET\k)+ 


(2.61) 




^gridJ Z r avg (k) 


i if v 0 ^ Zp • A/ • Oavgjk) 

1 + ? \v„ rid j 


H7\k-\)-Hp(k+{) 


2 \VgridJ 8 r,avg(k) 

Finally, the free-space electric field grid equation is 


EW)*E;-'m+Zo(j^) 


H’iHk-b-Hpi.k+h 


(2.62) 


Similarly, magnetic field update equations can be converted to magnetic field grid 
equations by replacing the variables z and t with the grid coordinates of the magnetic field 
spatial and temporal sampling points 
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z —» (&+^)A/, A: = 0, l...A^ z - 1 and t->(n + |)At, n = 0,1„ JV, - 1 
The magnetic field grid equation, using the same notation as for the electric field grid 

equation, is 

H,*kk+ ;) « - ) j—r-K(*+ !)-£;(*)] (2-63) 

or, for non-magnetic media 

Hf\k+\)»H7'{k+\)-y«y^][E",(k+l)-E’;m (2.64) 

The electric and magnetic field grid equations have the following general form 

E n r = C E \Ef + CnVHf (2.65) 

Hf w = C m Hf + CmVEf (2.66) 

where C's are constants (real numbers) that depend on the media properties and the 

velocity ratio v^v^, and "del" operator represents the spatial derivative (gradient). This 

general form of the grid equations can be interpreted as follows: "the new value of EJH 

field at a grid node is equal to the weighted sum of the old value of the E/H field at the 

same node and the spatial variation of the old H/E field between the two 

nearest-neighbor nodes ." 

2. Grid Termination 

The grid equations derived above are valid for all the nodes except for the nodes 
on the grid edges. The reason is that the grid edge nodes have only one neighbor node 
instead of two like any node that is not on the grid edge. The edge nodes with a single 
"neighbor" thus need to have equations different than the equations we have derived for 
the "non-edge" grid nodes with two "neighbors". TGT requirement is conceptually very 
simple and straightforward to implement in 1-D. To that effect we will use the concept of 
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a multiport (a two-port in 1-D). It does not matter, in principle, whether the edge node is 
an E field or an H field node. We will assume that the edge node for the two-port model 
is an E field node. The input port of the two-port will be the nearest node "of the same 
kind". Since we have assumed an E field node for the output port, the input port will be 
the nearest E field node inside the grid. We could have also selected an E field node as the 
output port and the nearest H-field node as the input port. However, the selection of the 
same kind of node for both input and output ports has the advantage that the TGT results 
obtained in this mann er can be also used to solve the wave equation (a second order 
partial differential equation) that has only one field as the variable and the grid with only 
one kind of nodes. The figure below shows the two ports for the two (1-D) grid edge 
nodes. 


out 


h(t) 


in 


-o 


o 


m 


o 


h(t> 


out 


grid edge 


grid edge 


Figure 2.6. Modeling 1-D Grid Terminations as Two-ports. 

The fields at the output ports E(0,t) and E(L,t) can be expressed as convolutions 
of the fields at the input ports E(Al,t) and E(L-Al,t) and the two port impulse response 
h(t). (The impulse responses for the two-ports on the left and on the right edges are 
identical, by symmetry). 

E x ( 0,0 = E x (Al, t) * h(t) = f 0 E x (Al, x)h(t - iyk (2.67) 

E X (L, t) = E X (L - A/, t) * h(t ) = f 0 E X (L - A/, x)h(t - x)dx (2.68) 
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The discretized forms of these equations, involving samples of the fields at t=nAt, are the 
discretized convolutions 

E n M = Z p p ZE p x {\)-h"~P (2.69) 

E”(Nz) = Z P P Z Ex(N z -1) • h n ~? (2.70) 

The convolution equations express the fields at the edges as weighted sums of the 

time histories of the fields "just inside" the grid. In that respect they are equations of the 

same type as the equations for the non-edge grid nodes, except that they involve, in 

general, summations with more terms. However, the impulse response h(t), as will be 

shown, is a rapidly converging (to zero) function which reduces the number of relevant 

terms in the convolution s um The impulse response h(t) (actually its "sampled" form h n ) 

needs to be determined only once, for a selected grid velocity v^^Al/At. The issue 

r emains how to determine the impulse response? The discretized impulse response h n will 

be determined using the discrete equivalent of the Dirac deha function which we will 

denote as S“. Since there are two grid edges and two input ports (Fig. 2.6) one input port 

will be set to zero and 8” will be applied to the other input port. Since, in 1-D, there is 

only one h n to determine we will set E x (Al,t) to zero and apply the Dirac delta function as 

E X (L-AU) 

E"(N Z - 1) = 8” (2.71) 

£?(!) = 0 (2.72) 

This is depicted in Figure 2.7. Note that the inpulse response h n is observed at z = 

L and that the grid extends, theoretically, to infinity past the observation point z = L. 

Since 8 n = 0 for n > 0, the grid to the left of the (L - Al) node is effectively isolated form 
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the observation point at L. This means that, to find h n , we may just consider the grid to 
the right of the L-Al node, as shown below. 

E( Al,t)=0 E(L-AU)=8(t) 

^ h(t) 

• O.•.O.•.O.•.O.•.O.•.o.•.O.•.o • 

A A 

gridec^e grideck*e 

E(L-Alt)= 5(t) 

•.e.•.o •.0.•o.• 

A 

gridafee 

Figure 2.7. Determining the Impulse Response. 

The impulse response needs to be obtained as if the grid were not terminated at all 
to the right of the observation point at z = L. If the grid were terminated, the 
"reflectionless" termination condition would be needed and this is exactly what we do not 
have and are trying to find. A grid extending to infinit y does not require a termination 
and thus "avoids" the termination problem. Although it is not possible, in practice, to 
extend the grid to infinity, the grid can be made large enough such that any "reflections" 
off the new grid edge would have arrived after the impulse response has "converged" to a 
small selected value that we consider as "zero". The time it takes the impulse to 

"propagate" to this new grid edge and back to the observation point is 

T d = ^- = ^^- = 2N D At ( 2 . 73 ) 

v grid A/ 

At 
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where D is the distance between the observation point and the new grid's edge and N D is 
the number of nodes between the observation point and the new grid's edge. The duration 
of the impulse response h“ that will not be "corrupted" by the reflections will be 2N D . 
Finally, it is assumed that the space outside the grid is ffee-space, and the grid equations 
for the electric and magnetic fields will thus be the free-space equations 

E;(k) = E;-\k) +Zo(^) H";\k-\)-H7'(k + \) (2.74) 

Hf\k + \)=Hf\k+\) + Yo(^][E%k)-Ei(k+l)\ (2.75) 

The order of the terms in the brackets for the magnetic field equation has been 

reversed, such that the electric and magnetic field grid equations will have the same sign 

(+) in front of the brackets. The above equations can be also written in a different form, 

to reduce the number of multiplication's that would need to be done at each time step. 

Multiplying the magnetic field equation by Z 0 we get 

ZoHy + hk+\) = ZoH n P(k+\) + (^) [E n x (k)-E"(k + 1)] (2.76) 

An auxiliaiy variable hy may be introduced 

h y = Z 0 H y (2.77) 

and the electric and magnetic field grid equations can be written using hy, instead of Hy, 

£;(*)= £r‘(*) + Ggj) h7'\k-^-h7\k + \) (2.78) 

kf\k + i) = hf‘(k + i) + \E%K)-E’(k + 1)1 (2.79) 

The equations for the electric field (E x ) and for the magnetic field multiplied by Z 0 
(Z 0 Hy) have identical forms, involving only one parameter, the velocity ratio v/v^. The 
boundary impulse response h n , to be obtained using these equations, will thus be valid only 
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for a selected velocity ratio or, equivalently for the selected grid velocity v^^Al/At The 
boundary impulse response h n , obtained using the equations for v ( /v grid =l, is shown below. 
Note that the plotting program interpolates linearly the values between the sampling points 
nAt. 



Time Step 

Figure 2.8. Boundary Impulse Response. 

The boundary impulse response for 1-D is very simple 

h n — S 1 (2.80) 

The impulse response is non-zero only for n = 1 where it has the value of 1, which is a 

delta impulse delayed by At and may also be written as 

h(t) = S(t-At) (2.81) 

We can now write the grid equations for 1-D grid edges as follows 

E”(0) = Z P p ZE p x (\)-h"-P =ET l (\)-h x =E n ~\ 1) (2.82) 

E n x (N z ) = 2^ E P X (N Z - 1) • h n ~ p = E n ~ x (N z - 1) • h' = E? l (N x - 1) (2.83) 
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The above equations state that the fields at the grid edges are updated by simply 
taking the previous values of their nearest neighbors inside the grid. Although the 
boundary impulse responses are not so simple for 2-D and 3-D, they can be determined 
using essentially the same procedure as shown for 1-D. 

C. 1-D TGT RESULTS 

We apply the geometrical model shown in Figure 2.6. The nodes on the edges are 
E field nodes. The source is the electric field at the center node and the source waveform 
is a unit amplitude delta pulse. This represents, in 1-D, a uniform plane wave propagating 
from the center into +z and -z directions. By applying the "standard" FD-TD equation for 
the nodes inside the grid and the TGT equations for the grid edges we obtain the power 
within the grid as a function of time as shown below. It is clear that the power within the 
grid is constant until the wave has "left" the grid when it falls to zero abruptly. The 1-D 
grid with TGT termination thus behaves like an "ideal" grid. 
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Figure 2.9. Residual Power for TGT. 

We will next compare this result with that for an "infinite" grid. Note that the 
power is calculated within the same grid region as in the previous figure. 



Figure 2.10. Residual Power for an " Infinit e" Grid. 
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m. ANALYSIS OF TGT FOR 2-D FDTD 





Figure 3.1. TE Z Field. 

The direction of propagation is indicated by the propagation velocity vector v . 
The unit vector in the direction of propagation and the electric and magnetic field vectors 

form a triplet of mutually orthogonal vectors. The "dual" TM z field has the components 

H*, Hy, and E z 

~H m(x,y, t) = H x (x,y, t)~x +H y (x,y, t)~y (3.4) 

E m(x,y, t) = E z (x,y, t)z (3.5) 

as shown below. 

H 

Hx .-.-.~> x 

Figure 3.2. TM Z Field. 
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Our objective is to determine the discretized forms of Maxwell's curl equations in 
the integral form for TE Z fields 

f Cg \_E x (x,y,t)~x +E y (x,y,t) y] • dl E = \i(x,y)H z (x,y, t)~z • ds E ■ (3.6) 

\_Hz(x,y, t)~z ] • dl H = zix^E^y, t)~x +E y (x,y, t)~y ] • ds H J+ (3.7) 

cj(x 5 y)[E x (x,y, t)~x +E y (x,y, t)~y ] • ds H 

and for TM Z fields 

j Cs \_Ez(x,y, t)~z ] • dl E = “| \i.(x,y^H x (x,y, t)x +H y (x,y,t)y ] • ds E ■ (3.8) 

f Cff \_H x {x,y,t)~x +H y (x,y,t)~y^* dl H = (3.9) 

e(x,y)[E z (x,y, t)~z ] • ds H ■ +JJ^ a(x,y)^E z (x,y, ffz ] *ds H 
Note that the the second curl equations do not have the source current terms, since 
it has been assumed that there were no source currents in the domain of interest. The 
contours of integration for the electric and the magnetic field circulations are, in general, 
different and are thus labeled C E for an electric field contour and C H for a magnetic field 
contour. Similarly, the surfaces associated with the contours are labeled S E for the 
magnetic flux and S H for the electric flux. The electric and magnetic field contours for a 
TE Z field are shown below. 
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Ex(x,y+dy/2,t) 

,. . 

Ey(x-dx/2,y,t ) i y _fEy(x t dx/2,y,t) 
*4iz(x,y,t)+. 

i ^ ij ; 

... 

Ex(x,y-dyj/2,t) 


Figure 3.3. Contours for TE Z Electric and Magnetic Field Circulations and Fluxes. 

The electric field contours C E and the magnetic field contours C H are square (A1 by 
Al) contours, but in orthogonal planes. The C E contours are in the xy-plane, while the C H 
contours are in the yz-plane. The contours can be compared to links of a (2-D) chain 
fence, evoking the idea of electric and magnetic field linkage in two orthogonal directions. 
The "dual" C E and C H contours for a TM Z field are shown below 



Hx(x,y-dy/2,t) 


Figure 3.4. Contours for TM Z Electric and Magnetic Field Circulations and Fluxes. 




























The electric and magnetic contours and their associated surfaces will be used to discretize 
Maxwell's curl equations. 


2. Derivation of TE Z Magnetic Field Update Equation 

The discrete equivalent of the TE Z electric field circulation will he determined next. 
A contour C E is shown below 


Ex(x,y+dy/2,t) 


Ey(x-dx/2,y,t) 


JfS I Ey(x+dx/2,y,t) 

I Hz(x,y,t) I 



Ex(x,y-dy/2,t) 


Figure 3.5. A Contour for TE Z Electric Field Circulation and Magnetic Flux Calculations. 

The center of the surface S E is assumed to have the coordinates (x, y). We will 
assume a unif orm grid with Ax=Ay=Al. A local coordinate system (£, vj/) will be 
estabhshed, with the origin at the center of the contour such that any point (x',y*) within or 

on the contour C E can be specified by its local coordinates \ and \\i 

x' =x + £, y'=y + \y where - y < £ < +^- and -~<\y< +y 

The local coordinates will be used in evaluation of the line and surface integrals 

that constitute the integral forms of Maxwell's equations. The circulation of the electric 

field around the electric field contour C E , assuming counter-clockwise reference direction 

-4 

such that the normal to the surface S E is in the direction of the magnetic field H(x,y 7 i ), 
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can be evaluated approximately, assuming that the electric fields are constant over the 
length A1, along the edges of the (A1 by Al) square contour 

j Ce \E x (x'y, f)~x +Ey(x',/, t)7 ] • dl E « (3.10) 

[E x (x,y - f, 0 +E y (x + f,y, t) - E x (x,yE y (x - f,y, f)]A/ 

The magnetic flux through the surface S E can be calculated using the local coordinates 

JJ & ()?•*£= (3.11) 

G 0 (ifr+5.>'+ +%,y + 'V, ')"? • 

There are infinit ely many ways to model the variation of the magnetic flux density 

with the local coordinates \ and vj; within the contour. We first need to postulate a 

certain variation of the magnetic field with | and vp such that the integral can be evaluated 

over S E . The simplest model assumes that the contour width Al is small enough such that 

the magnetic field H(x+£,y+i|/,t) within the contour may be assumed constant and equal to 

the value at the contour's center H(x,y,t). 

H z (x+^,y+\\f,t)»H z (x,y,t) (3.12) 

This is equivalent to using a piece-wise constant or 2-D "pulse" expansion 

approximation of the actual magnetic field variation with the z-coordinate. The above 

assumption allows that magnetic field , although constant within a contour, can change 

from one contour to an other. This yields an approximate expression for the for the 

magnetic flux 

ff \l(x' y)H z {x r ,/, t)~z • ds E * Hz(x,y, t) • J + | 0 p(x + $,y + y)d^dy (3.13) 

JJ Sg 2 2 

The integral of the permeability p(x+^,y+vp) can be re-written in the following manner 
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0 0 n(x + %,y+vWv = (AO 2 • —0 0 0 v(x++ vy&iy (3.14) 

2 2 |_ (A/) 2 2 

The term in the brackets is recognized as the average permeability within the contour 

C E 

= 7-72 O 0 ►*(* + ^ + < 315 ) 
(A l) 2 2 

The approximate expression for the magnetic flux through S E can now be written using the 
average permeability as 

|0 \x(x' ,y')H z (x' ,y', t)~z • ds E ~ (M) 2 \iavg(x,y)H z (x,y, t) (3.16) 

Again, this approximate expression resulted from the piece-wise constant 
approximation of the magnetic field with respect to x and y-coordinates. The main 
advantage of the piece-wise constant expansion employed above is its simplicity. Better 
accuracy can be achieved by using more involved models for the field variation with x and 
y, but at the expense of increasing the computational time. The first curl equation of 

Maxwell can now be replaced by an approximate equation 

\_E x (x,y - f,0 + E y (x + f,y,i) -E x (x,y + E y (x - f,y, o] • A/« (3.17) 

-§ { [V«vg(x,y)H z (x,y, OKA/) 2 } 

which can be simplified as below because of media stationarity 

E x (x,y-f,t)+E y (x + f,y,t)-E x (x,y + f,t)-E y (x-f,y,t)* (3.18) 

-J \{Hz{x,y, 0) • A/ • Pov g (x,y)} 

This equation can also be re-written as 

|{// z (x,y,0} (3-19) 

A/-p~|.(y y) [^( X;> ’“ ^ +Ey(?C + **~ Ex(x ’ y Ey(x ~%’ y ’*)] 
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The time derivative operator in the above equation is typically "replaced" by the 
finit e difference approximation. However, just like in 1-D, we present an alternate 
approach, such that the approximation of the field time variation is shown to be analogous 
to the approximations already introduced for the field spatial variation. The above 
equation is integrated with respect to time to get an approximate equation for the 


magnetic field at the present time t 




(3.20) 


[ £ E,(x,y - f, i )* - £ E,(x,y +f ,*)*+!'„ E,(.x + f-H *>* - £ E A X ~ f T )*] 

A similar integral equation can be written for the magnetic field at the same spatial 


location but at an earlier time t-At 




(3.21) 


E x (x,y - f, x)dx - E,(x,y+%, x)dx + ' E y (x + f,y, x)dx - £* E y (x - f,y, i)dx] 


Subtracting the two equations we get 


H z (x,y, t) * H s (x,y, t-At) - 


(3.22) 


Al-Hav&y) 

[^E,(x,y-f, Z yh-^E^,y+f,z)ch+f t _ u E,(x+f,y,Tyk-f t _ it E,(x-f,y,x)dz] 
The integrals on the right-hand side can not be evaluated exactly, because the 

exact temporal variation of the electric fields within the At interval prior to t is generally 

not known. However, the integrals can be evaluated approximate^ by assuming a certain 


variation of the electric field with the temporal variable t. The amplest approach, 


consistent with the assumptions made for the field spatial variation, would be to assume 
that the interval At is small enough such that the electric field can be assumed 
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approximately constant within At and equal to the value at the center of the time interval 
(t-At, t) 


E*(x,y - f > T > * E &<y ~ f ><- f) 

for 

t-At <T < t 

(3.23) 

E,(x,y+f,x)*E,(x,y + f,t-%) 

for 

t- At <x <t 

(3.24) 

3k* 

I 

PS 

1 

II 

e 

PS 

1 

for 

t- At<T <t 

(3.25) 

E y (x+ §,y, 7) **Ey(x+f,y,t- f) 

for 

t-At<x<t 

(3.26) 


The above represents a piece-wise constant approximation of the electric field variation 
with respect to the temporal variable t. The approximate expression for the magnetic field 


at location (x,y) and at the time t now becomes 

(3.27) 

The above equation can be written in a somewhat different form, using the identity 

_1_1_ Jzo_ _1_1_ 1 _ 1 1 _1_L _ v y n J_ (2 281 

APSES’ fiAjri^ ZoV-r P-r 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Z 0 « 377H), and p r denotes relative permittivity. Introducing the 
grid "propagation" velocity 


_ A/ 

V grid ~ Af 


(3.29) 


we can write the approximate "update" equation for the magnetic field as 
H,(x,y,0 * H,(x,y,t-At )-x 


(3.30) 


[E x (x,y-f,t-j)- E x (x,y + f,t-f)+E y (x + f,y, t-f)-E y (x- f ,y, t- f ] 


The equation simplifies for the case of non-magnetic media (p f — 1) 


H z (x,y, t) « H z (x,y, t-At) - To^ x 


(3.31) 


[£,(*,>- f / - f) - E„(x,y + f , t - f) +E y (x + f ,y, I- f) -E y (x -f,y,t-f)] 
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The last two equations show that the magnetic field H and the electric field E are 
evaluated at points shifted spatially by Al/2, and at instants separated in time by At/2, just 
like in 1-D. The relationship between spatial and temporal "samples" of the electric and 
the magnetic fields is thus the same and the samples are shifted with respect to each other 
by one-half of the sampling interval. The approximate equations for the electric field 
"updates" can be derived using the same procedure shown for the magnetic field update 
equations. 


3. Derivation of TE Z Electric Field Update Equations 

We will start with the Maxwell's curl equation for the circulation of the magnetic 


field 


H z (x,y, t)~z • dl H = J*; j jJ Sw z( x ,y)\_E*( x ,y, if* +£y(x,y, 0 y ] • ^sh 
JJ Sw o(x,y) E x (x,y,t)x +E y (x,y, fyy ] • ds H 


+ (3.32) 


The discrete equivalent of the magnetic field circulation needs to be determined 
for two sample contours a contour in a plane parallel to the yz-coordinate plane, and a 
contour in a plane parallel to xz-coordinate plane. 



Hz(x-A l/2,y,t) 


0/ “ Ey(x.y.,)^ 


> * 


Hz(x + Al/2,y,t) 


Figure 3.6. A Contour for Updating E y . 
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A local coordinate system (£, Q will be used, with the origin at the center of the 
contour. Any point (x',z') within or on the contour C H can be specified by its local 


coordinates 


x A =x + £ z' = z+C, where - — <£,< +— and 


2 ~^~^2 


The circulation of the magnetic field around C H , assuming counter-clockwise reference 

direction such that the normal to the surface S H is in the direction of the electric field 
E y (x,y, t), is therefore given exactly by the following simple expression (because there is 
no magnetic field variation with the z-coordinate) 

f Hzix'y, ijt *dh = 0 H Z (X - f ,y, f)~z • d£z - 0 H : {x + f ,y, i)z • d^z (3.33) 


0 H z (x',y', t)z.dl H = [H z (x - f ,y, t) - H z (x + f ,y, f)] - A/ (3.34) 

The electric flux through the surface S H can only be evaluated approximately, since 
the exact way that the electric flux density D y (x,y, t) = z(x,y)E y (x,y, t) varies with x and y 
is not known a priori. First we need to evaluate the flux integrals 


J0 K x ',y'd) E y(x',y' , t)~y • ds H = 
0 0 e(*+ z>,y ++ £»y +v, 0 ~y • 


(3.35) 


0 o(x , ,y')Ey(x'y, t)~y •ds H = 

0 0 c ( x+ + v) E *( x + z»y +v, t)~y • ~y&& 


3.36) 


Since the integrands are not functions of the local z-coordinate C, and the local 


coordinate vj/=0, the surface integrals are simplified to line integrals 

J0 s(x+ \,y)E y (x + %,y, t)~y •ds H = Al 0 e(x+ ,y)E y (x + \,y, t)d£, (3.37) 
IT cr(x + l,y)E y {x + l,y, t)~y • ds H = A/ 0 a(x + B,,y)E y (x + B^y, t)dE, (3.38) 
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Next, we assume that the contour side A1 is small enough such that the electric 
field E(x+£,y,t) within the contour may be assumed constant and equal to the value at the 
contour's center E(x,y,t). This assumption yields approximate expressions for the flux 


integrals 

[T s(x + %,y)E y (x+%,y, t)~y • ~ds H ~ A/ • E y (x,y, t) • f + J s(x + ^,y)d^ (3.39) 

f|" a(x + l,y)E y (x + Z s ,y,t)y • ds H * A/• E y (x,y,t) -0 a(x+ ^y)dt (3.40) 

JJ S H 2 

or, using the permittivity and conductivity averaged in the x-direction 

IT s(x + $,y)E y (x + l,y,i)y • ds„ * (A/) 2 -E^yjy-e^iyy) (3.41) 


IT o(x + y)Ey(pc + j 7 , t) y • ds h « (A/) -E y (x^y^ t) * <5avgx(x,y) P-42) 

JJ s H 

The second curl equation of Maxwell can now be replaced by an approximate equation 

[H z (x - f,y, t) - H z (x + f ,y, t)] ■ A / * (3.43) 

j t {[s a Vg x(x,y)E y (x,y, t)](M ) 2 } + (A if • <w(x,y) • E y (x,y, t) 

which can be simplified to 

H z (x - f ,y, t) - H z (x+ f ,y, t) « (3.44) 

£ avgx(x. _}) • A/ • {E y (x,y, t) ) + A/ • X5avgx(y r y) ■ E y (x,y,f) 


The approximate equation can now be re-written as 

~ (3.45) 

.. *, - f,y, 0 - H,(x + f,y, 4)] - • £,(*,* i) 

■Zavgx(x,y) L 2 Zavgx(X,y) 

The above equation is integrated with respect to time to get an approximate equation for 


■E y (x,y,t) 


the electric field at the present time t 
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E y (x,y, t)« — Sgv ^ (x ;^ [f 0 H ^ x ~ H z (x + f,y,x)di] (3.46) 

c Javgx(x,y) f r £-i /_ _ 

A similar equation can be written for the electric field at the same spatial location but at an 
earlier time t-At 

E y (x,y, t- AO * A/ , £ y r H &- f*y > T )^ - r ^ z(x + f > ^ T ^ T ] ( 3 - 47 ) 

_ 

Wfcj) J o * v 


Subtracting the two equations we get 

E y (x,y, t) -E y (x,y,t-At) * 

■KTI^yl 

Gavgxfay) F p , w 


(3.48) 


which can also be written as 


E y (x,y,t) 


(3.49) 


'" ■ AO '+ A/.g^j L H * X ~ f’ y ’ X)<h ~ Lr ^ z(x + * T H 

-^4f Ey(x,y,x)ch 
Savgtfry) ^ 

If the interval At is small, the magnetic fields can be assumed approximately constant 


within At and equal to the value at the center of the time interval At 


H z (x+f,y,T)*H z (x + f,y,t-f) for t-At<x<t 
H,(x-%,y,t)*H.(x-%,y,t-f) for r-ArStir 


(3.50) 

(3.51) 


Using the above approximation we get 
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(3.52) 


E y (x,y,t-At) + 


E y (x,y, t) » 

A ' {ft<* - f ,y, t - f) - H,(x +f t - f)] 


A/ • Savgx(x,y) 


^avgx^y) f r p, 

The above equation can be written in a somewhat different form 

E y (x,y, t )* 


(3.53) 


E y (x,y,t- At) + 


At 


A l • Savgxfay) 


[H,(x - f ,y, t - f) - ft<* + f t - f)] 


The term in the brackets is recognized as the average value of the electric field component 

E y at (x,y) within the interval (t-At, t). The above equation can be written in a more 

compact form using the identity 
1 1 1 


1 


1 


1 


_____ ft! _1_ 7 J_ 

s So Sr jm^VZotr °°Sr 


(3.54) 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Zo « 3770) and s r denotes relative permittivity. Introducing the grid 
"propagation" velocity 

< 3 «) 

we can write the approximate electric field "update" equation as 

E y {x,y, t) » (3.56) 

E y (x,y, t- At) + Zo v ° rid Zrav ^ X y) [Pz(x - f ,y, t- f) - H z (x + f ,y, t - f)] 

" Z °vJride r *J'(x,y) * A/ ' ^^^y)Ey^)(x,y, t) 

The average of the electric field between t-At and t, denoted E x>Jvg(At) can not be 


calculated exactly because the exact temporal variation of the electric field within the At 
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interval is not known a priori. However, assuming a linear variation within At interval, the 


average field is the average of the values at t-At and t 


Ey,avg(&t) (x, 0 1 


E y (x,y, t- At) +E y (x,y, t) 


(3.57) 


Substituting the time-average E y avg(At) (x,y,t) into the update equation we get 


E y (x,y, t) « (3.58) 

Ey(x,y, t-At)+Z 0 v g l dSravg ](x,y)[ Hz ( X ~ %’ y ’ H *( x + * y > f ~ 

^ v 0 1 A , _ Ey(x,y,t-At)+E y (x,y,t) 

z °^Z^M' Al ' Gavgx(x ' y) 2 

which gives the final equation for the electric field updates as 

1 1 Vp Zq • A/ • g a vgx(x, j) 

E y (x,y, t >« 2VPU *'Ta' y ]xv^ ’ y - ‘" A<)+ (3 ' 59> 

l 1 v 0 ^0 * A/ • Oqvgx\X,y) 

2 V grid £r,avgX(x^y) 

Z v ° 1 

Vgrid &r,avgx(x,y) 

r TT / AZ f TT ( . tJ f AK~| 

- —————-h # z (x - t) ~ h a x + ~z>y> * ~ t'J 

1 | 1 Vo Zp-Al- Gavgx&y) L 

2 Vgr’d £ raV gx(x,y) 

Hie equation simplifies greatly for the case of non-conductive media (cr=0) 

E y (x,y, t) « (3.60) 

Ey(*,y, t - At)+Zo Vg ° rid Sr av ^ y )[ H -( x ~ f *y> 1 - f) - H *( x+ f >y> *-t>] 

In the case of free-space (e =1), the equation simplifies further 

E,(x,y, I ) * E y (x,y, I - Ar) +Z 0 *j*j[A(r - f ,y, / - f) -ft(x+fO-, f- f)] (3.61) 

Hie above equations are the update equations for the y-directed component of a TE Z 

electric field. S imilar equations will be formulated next for the x-directed electric field 


(3.60) 


component. 
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A 




Hzj(x,y+dl/2,t) 


j Ex(x,y,t) 




Hz(x,y-dl/2,t) 


Figure 3.7. A Contour for Updating E x . 

A local coordinate system will be used, with the origin at the center of the 
contour. Any point (y'.z 1 ) within or on the contour C H can be specified by its local 
coordinates 

y' =y+ \y z' = z+Q where -^<vj/<+^ and 
Then, by applying the same procedure drived for E y field component, we get the final 
equation for E x electric field component as shown below 


J 1 Vp Zq • A/ • Paygyipc^y) 

j-i , A 2 V g" d Zrj3vgy(x,y) F . . 

EJx.y, t) «--—77-7— r^x(x,yj - i\T)+ 

, J Vq Zq • A/ • CTtrygyC^y) 

+ 2*W Zr fi vgy(x,y) 

2 q v o 1 

-+f,f) - H,(x,y - f / - f)] 

1 Vo 2/p • A/ • CJavgy(X,y) 

2 £r,avgy(x,y) 

The equation simplifies greatly for the case of non-conductive media (cr=0) 


(3.62) 
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(3.63) 


E x (x,y, t )» 

E x (x,y, t-At)+ Z ^v]° ridSrav ^(x,y)^ Hz(x, y + " f > " ~ aM '~ f>] 

In the case of free-space (s =1), the equation simplifies further 

E x (x,y, t) « E x (x,y, t - At)+Z 0 ^_H z (x,y + f, t - f) - H z (x,y - f,t- f)] (3.64) 

The above equations are the update equations for the x-directed component of a TE Z 
electric field. 

4. Derivation of TM Z Electric Field Update Equation 

The discrete equivalent of the TM Z magnetic field circulation will be determined 
next. A sample contour C H that will be used to determine the electric field update 
equation is shown below 


Hx(x,y+dy/2,t) 


Hy(x-dx/2,y,t) 


Ez(x,y,t) 


Hy(x+dx/2,y,t) 


Hx(x,y-dy/2,t) 

Figure 3.8. A Contour for TM Z Magnetic Field Circulation and Electric Flux 

Calculations. 

The center of the surface S H is assumed to have the coordinates (x, y). We will 
assume a uniform grid with Ax=Ay=Al. A local coordinate system (£,, vp) will be 
established, with the origin at the center of the contour. Any point (x , ,y I ) within or on the 

contour C E can be specified by its local coordinates \ and \\i 
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x / = x + £ /=j + ip where and ~Y<\j/<+— 

The local coordinates will be used in evaluation of the line and surface integrals that 

constitute the integral forms of Maxwell's equations. The circulation of the magnetic field 

around the magnetic field contour C H , assuming counter-clockwise reference direction 

—^ 

such that the normal to the surface S E is in the direction of the electric field E (x,y, f), can 
be evaluated approximately, assuming that the magnetic fields are constant over the length 

A1, along the edges of the (A1 by Al) square contour 

f Cw [#,(*',/, t)x +H y (x'y ,« (3-65) 

[H x (x,y - f, f)+ H y (x + f,y, t) - H x (x,y + f, t)- H y (x - f,y, f)] Al 
The electric flux (the "displacement current") through the surface S H and the 

current flux (the conduction current) can be calculated using the local coordinates 

JT e(x',y , )E z (x',y', t)~z • ds H = (3.66) 

fl G & ( x + + v) E z( x +%»y + o? * 

2 J 2 

IT o(x / ,y / )E z (x / ,y / , t)~z • dsn = (3 67) 

JJSff 

G G a ( x+ ^y + v) E *( x + *»y +v, 

2 J 2 

There are infinitely many ways to model the variation of the electric flux density 

with the local coordinates E, and vp within the contour. We first need to postulate a certain 

variation of the electric field with ^ and vp such that the integral can be evaluated over S H . 

The simplest model assumes that the contour width Al is small enough such that the 

electric field E(x+£,y+vp,t) within the contour may be assumed constant and equal to the 

value at the contour's center E(x,y,t) 

E z (x + \,y +ip, i) » E z (x,y, t) (3.68) 
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This is equivalent to using a piece-wise constant approximation of the actual electric field 
variation with the z-coordinate.The above assumption allows that, although the electric 
field is constant within a contour, it can change from one contour to an other. This yields 
approximate expressions for displacement and conduction currents 

ff 8 (x / ,/)£ r (x / y, ffz *~ds E « E z (x,y, t) • 0 0 e(* + %,y + M(3-69) 

f>H 2 2 

ff <j(x',/) E z (x',y', ffz • ~ds E « E z (x,y,t) • 0 0 a(x + + y)^U/ (3.70) 

The integrals of the permittivity s(x+£,y+vp) and the conductivity cr(x+£,y+y) can be 
re-written in the following manner 

0 0 e(x + i,y + y)^¥ = (AO 2 • —0 0 0 s(x + 1,7 + V)^V (3.71) 

T 2 L(A/) 2 2 

0 0 a(x + £,y + v|/)c^^vj/ = (A if ■ 7777 0 0 ct(x +^,y + vj i)dt,d\y (3.72) 

J ~2 2 L(A l) 2 2 

The term in the brackets is recognized as the average permittivity and the 
average conductivity within the contour C H . 

Stfv g (x,T) = 77-7 0 0 z(x+Z,,y+ ( 3 . 73 ) 

(A/) 2 2 

<Wx,>>) = 7777 0 0 cr (x+^,y +(3.74) 

(A/) 2 2 

The approximate expression for the displacement and conduction currents through 

S H can now be written using the average permittivity and the average conductivity as 

ff 8 (x' ,y')E z (x' ,y', ffz • ds H « (Al) 2 s avg (x,y)E z (x,y, i ) (3.75) 

JJ s H 

ff a{x'y)E z (x'y, iff • ds H « (A/) 2 a a vg(x,y)£ r (x,y, t) (3.76) 
JJ s H 

The first cirri equation of Maxwell can now he replaced by an approximate equation 
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[H x (x,y-fj)+H y <x+f,y,t)-H x (x,y+f,Q-H y (x-f,y,t)]AI* (3.77) 
i- l [lf.c,.' S (^.y)E ! (x.y, t)](Al) 2 } + [c,. t (x,y)E,_(x,y, 


winch can be Amp lified (because of media stationarity) to 

H x (x,y - f, t) +H y (x + f,y, t) - H x (x,y + H y (x - f ,y, t) * (3.78) 

| ] -{E z (x,y, /)} • A/ • Sav g (x,y) + a„v g (x,y) • E z (x,y, t) • A/ 

TM> equation can also be re-written as 

^ t {E,(x,y,t)}« (3.79) 

—— — — -[ff,(x,y -j,t)+H,(x+ §,y, I) - H,(x,y +f, t) - H,(x - j,y, f)] 
A/-e^(x,y) L 

&avg(x,y) p , a 

-7^y) E ^ 

The above equation is integrated with respect to time to get an approximate equation for 


(3.79) 


the electric field at the present time t 


*#•*’*> M u-'ljM* <3 80) 

[f o H,(x,y - f, iyk - J' HAx,y + f, iyk + {' H y (x +f.v, !)*-{' H,(x - f,y, x)ih\ 

-^jy^ EAw)dl 

A gimilar integral equation can be written for the electric field at the same spatial 
location but at an earlier time t-At 

E z (x,y,t- At) « (3.81) 

Hx(X ' y V HAx ’ y + T)A + C + 2 T >*- 

r ^ - i.y. t)*l - • c*>* 

Subtracting the two equations we get 
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(3.82) 


Ez(x,y, t) « E z (x,y, t-At) + 


1 


Al ■ Zavg(X,y) 


[L, H x (x,y - f, x)dx -H x (x,y+ f, x)dx +{^ H y (x + f,y, x)dx - tf,(x - f,y, z)dx] 




The integrals can be evaluated approximately by assuming a certain variation of 


magnetic field with the temporal variable x. The simplest approach, consistent with the 
assumptions made for the field spatial variation, would be to assume that the interval At is 
email enough such that the magnetic field can be assumed approximately constant within 
At and equal to the value at the center of the time interval (t-At, t) 


HAx.y - f. t)« H,(x,y 

for 

t-At<x<t 

(3.83) 

H,(x,y+f,x)*H,(x,y+f,l-f) 

for 

t-At<x<t 

(3.84) 

H y <x-f,y,x)«H y (x-f,y,t-f) 

for 

t-At<x<t 

(3.85) 

H y (x + f,y,t)*H y (x + f,y,l-f) 

for 

t-At< x < t 

(3.86) 


The integral of the electric field is recognized as the average value of the field 


within the At interval 


S avg(x,y) 




Povgfoy) • At 
s avg (x,y) 


' Ez,avg(6t)(x,y, t) 


(3.87) 


The average value of the electric field is approximately 


Ez,avg(At)(x,y, t) '■ 


E z {x,y, t) +E z (x,y, t-At) 


(3.88) 


The approximate expression for the electric field at location (x,y) and at the time t 


now becomes 


E z (x,y,t) * E z (x,y,t-At) + 


At 


Al * &avg(x,y) 


(3.89) 


[H,(x,y - f, t - f) - + f, t-- f) +H y (x + f ,y,t- f) ■- HJx - f ,y, 1- f)] 


Gcn,g(x,y) • A t E z (x,y, t) + E z (x,y, t- At) 

Zav g (x,y) ’ 2 
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The above equation can be written in a somewhat different form, using the identity 

i n = j^ 1 1 _ i [E* j_ = v z j_ 

5 So 5 '- 0 °Er 


(3.90) 

where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Z 0 « 377fi), and s r denotes relative permittivity. Introducing the grid 
"propagation" velocity 

= (3.91) 

we can write the electric field equation more compactly 

El (x,y,t) . E,(x,y,I-A<) + ^^^35 x (3.92) 


_ Vo Oavg&y) • A/ E z (x,y, t) + E t (x,y, t - A t) 

^ grid £ raV g (x. J/) 2 

The electric field "update" equation is obtained by grouping the like terms in the 


above equation 


1 1 Vq Zq ' Gavgfay) * 

~ 2Vgnd S navg&y) 

^ ^ 1 Vq Zo • Gavg( x,y) * A/ 


E z (x,y,t) 


E z (x,y,t-At) + 


(3.93) 


Zo 


Vo 


1 


V grid£ navg (x,y) 


. I Vo Zo-Oavgix^-Al 

^ I A 1 


' 2 V grid S rw&y) * ' 2 V W 8 r>a vg(x,y) 

•[ft(x,y - f, / - f) - H,(x,y + f, / - f) +f, * 1 ■- f) ■ -H,(x- *y,it -- f)] 

The equation simp lifies for the case of non-conductive media (c=0) 


vo 


E z (x,y, t ) a E z (x,y, t-At) + Z 0 — 


(3.94) 


V gridZ ravg (x,y) 

[H x (x,y - f , t - f ) - H x (x,y + f , t-- f )+ H y (x+ f ,y, t- f ) - - f ,y, i - f)] 

The approximate "update" equations for TMz magnetic field components will be derived 


next. 
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5. Derivation of TM i Magnetic Field Update Equations 

We will start with the Maxwell's curl equation for the circulation of the electric 


E,(x,y , ijt • dig- -J~ s(x,y)^H x (x,y, f)~x +H y (x,y, t)y ']•<**■ (3.95) 
The discrete equivalent of the magnetic field circulation needs to be determined 
for two sample contours a contour in a plane parallel to the yz-coordinate plane, and a 
contour in a plane parallel to xz-coordinate plane. 


i t 

Ez(x-Al/2,y,t^/ ^ Hy(x,y,t)^/^ 


Ez(x + Al/2,y,t) 


Figure 3.9. A Contour for Updating Hy. 

A local coordinate system (£, Q will be used, with the origin at the center of the 
contour. Any point (x^z 1 ) within or on the contour C E can be specified by its local 

coordinates 

x' =x + 't > z' = z+C, where ~^~ <^< and -^<£<-1 — 

z* z* z* z* 

The circulation of electric field around C E , assuming counter-clockwise reference 


direction such that the normal to the surface S E is in the direction of magnetic field 
H y (x,y, t), is therefore given exactly by the following simple expression because there is 
no electric field variation with the z-coordinate 
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( 3 . 96 ) 


j Cg E z (x',y', t)~z • dl E = 

fa E ^ x -f>y> • d &- 0 Ez(x + f >y>Q* • 

2 2 

j Cg Etix'y,(ft •alt- [E z (x -f ,y,t)-E z (x. + f ,y ,/)] • A/ (3.97) 

The magnetic flux through the surface S E can only be evaluated approximately, 

since the exact way that the magnetic flux density 

B y (x,y, t) = p(x,y)H y (x,y, t) (3.98) 

varies with x and y is not known a priori. First we need to evaluate the flux integral 

JJ Sg I x(x',y', t)H y (x',y', t)~y • ds E = (3-99) 

0 0 ^( x + ^y + v) H y( x+ z»y + v, t)y • ~y 
2 2 

Since the integrand is not a function of the local z-coordinate £ and the local coordinate 
vj/=0, the surface integral simplifies to a line integral 

IT p(x + l,y)H y (x + ^,y,t)y »ds E = A/ fa [i(x + ,y)H y (x + E,,y, t)d^ (3.100) 

Se 2 

Next, ass ume that the contour side A1 is small enough such that the magnetic field 
H(x+£,y,t) within the contour may be assumed constant and equal to the value at the 
contour's center H(x,y,t). This assumption yields an approximate expression for the flux 
integral 

IT |i(x + ^,y)Hy(x+t)~y • ds E ~ A/- H y (x,y,t) • fa p(x + £,y)d£ (3.101) 

Se 2 

or, using permeability averaged in the x-direction 

fa p(x + l,y)H y (x + t)~y • ~ds E * (A/) 2 • H y {x,y, t) • p^(x,y) (3.102) 

Se 

The first curl equation of Maxwell can now be replaced by an approximate equation 

[£,(x - f ,y, /) - E-.(x + f,y, <)] • 4/ » -f,{l^(x,y)H 1 ,(x,y, /)](A/) 2 } (3.103) 
which can be simp lified to below equation, because of media stationarity 
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E z (x - f ,y, t) - E z (x + f,y, t) * ~[i mgx (x,y) ■ A / • E-{H y (x,y, t ) (3.104) 


The approximate equation can be re-written as 

|{^(x,y, /» * jr -±^ [E 1 (x - f■ * 0 - Ei(x + fo', 0] (3-105) 

The above equation is integrated with respect to time to get an approximate 
equation for the electric field at the present time t 

Hy(x,y, t) * A/ ,^; (y -^[lo Ez(x ~ fj’ x)dx ~ Jo Ez(x + V' y ’ x)dx ] (3 ‘ 106) 

A similar equation can be written for the electric field at the same spatial location 
but at an earlier time t-At 

H y (x,y,t-At)* A/ .^ [C E *( x ~ f’T,-£* E z (x + f ,y,x)dx\ (3.107) 

Subtracting the two equations we get 

H y (x,y, t) - H y (x,y, t-At)* (3.108) 


A/-Hovs*(x,t) 
which can also be written as 


■ ullfaJ L, ~ ** T) * - L + f-y- ')*_ 


H y (pc,y, t) * (3.109) 

H r (x,y,, - A0 + E z (x-%,y,xyk-E z (x + f,y, t)*] 

Assuming that the interval At is small enough such that the electric fields can be 

assumed approximately constant within At and equal to the value at the center of the time 


interval At 


E z (x + f,y,x)*E z (x+f,y,t- f) for t-At<x<t 

£J * £J\ 


we get 


E z (x- T ,y,x)*E z (x- T ,y,t- T ) for t-At<x<t 


H y {x,y,t)* 

H„(x,y, t-At) + + f,y, t - f) - E,(x - f ,7, < - f)] 


(3.110) 

(3.111) 

(3.112) 
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The term in the brackets is recognized as the average value of the electric field 
component E y at (x,y) within the interval (t-At, t). The above equation can be written in a 
more compact form usingthe identity 

1 _ 1 1 _ _1_1_ _ 1 / So 1 _ v y 1 /-j in\ 

V~V0llr~ ° ° Hr V ' ' 

where the free space velocity of propagation is denoted v 0 , the intrinsic admitance of free 

space is denoted Y 0 (To = -~~,Zo «377Q) and p r denotes relative permeabihty. 

Zo 

Introducing the grid "propagation" velocity 

(3.114) 

we can write the approximate equation for the electric field "update" as 

H y (x,y,t)* (3.115) 

In the case of free-space (|o.=l), the equation simplifies to 

H y (x,y, t) w (3.116) 

H y (x,y, / - At) + + f ,y, t-f)-E,(x- f ,y, t - f )] 

The above equations are the update equations for the y-directed component of a TM Z 

magnetic field. S imilar equations will be formulated next for the x-directed magnetic field 

component. 
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Ez(x,y+dl/2,t) 


—Hx(x,y,t) 




Ez(x,y-dl/2,t) 


Figure 3.10. A Contour for Updating 


A local coordinate system will be used, with the origin at the center of the 
contour. Any point (y*,^) within or on the contour C E can be specified by its local 


coordinates 


y'=y + \\) z' =z+C, where -^<v|/<+^ and <t,< +^j- 

By following the same procedure shown for Hy update equation, we can determine 

the approximate equation for the H* magnetic field "updates" as 

H x (x,y, t) « (3.117) 

H x (x,y, t-At)+ y 0 { E &*y -f’t-f)- E *( x >y + ■~ f) } 

In the case of free-space (|i=l), the equation simplifies to 

H x (x,y,i)^H,{x,y,t-M) + Yo^[E z (x,y-^,l-f,-E z {x,y+^t-^)\ (3 118) 
The above equations are the update equations for the x-directed component of a TE Z 

electric field. 
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B. TRANSPARENT GRID TERMINATION IN 2-D 


1. 2-D TE I and TM I Grids 

Hie electric and magnetic field update equations have been derived using local 
coordinate systems, with ori gins at the centers of the magnetic and electric contours, 
respectively. These equations now need to be "converted" to a global coordinate system, 
that is to the grid of equi-distant sampling points in the xy-plane. Since we have either 
TE Z or TM Z 2-D electromagnetic fields there will be TE Z and TM Z grids as well. Let us 
assume that our do main is an L by L square. We will assume a uniform grid, that is a 
discretization step A1 that is not changing with position. The fields are then "sampled" 
using a spatial step A1=L/(N-1). The electric and magnetic field sample locations are 
"interleaved", as shown below for N=5. The spatial grid edge samples can be, in principle, 
either electric field samples or magnetic field samples. We will select electric fields for 
spatial edge samples for TE Z grids and magnetic fields for TM Z grids. The "generic" 2-D 
grid below is thus applicable to either TE Z or TM Z fields. 
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Figure 3.11.2-D Grid 

An (N by N) TE Z grid has (N)(N) magnetic field nodes, (N-2)(N-1) E x electric 
field nodes, and (N-l)(N-2) E y electric field nodes. Similarly, an (N by N) TM Z grid has 
(N-2)(N-1) E z electric field nodes, (N-l)(N-2) H z magnetic field nodes, and (N+1)(N) 
magnetic field nodes. TE Z electric and magnetic field "update” equations can now be 
"converted” to electric and magnetic field grid "update" equations by replacing the 
variables x, y and t with the grid coordinates of the electric field spatial and temporal 
sampling points 

x ->■ /A/, i = 0, 1...N y -> jAl, i = 0, 1...N t -> nAt, n = 0,1..JV, 
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The notation can be simplified by omitting the common A1 and At terms, and using 


a superscript for the index of the temporal sampling point. The grid equations for 


conductive media are obtained from equations (3.62,3.65) 

j 1 Vo Zo ' A/ • (T avgyijij) 

m,j)m \^ iJ) -sr'M 

J 1 Vo * A/ * CT avgyijij) 

2 Vgrid £ r,avgy(J)j) 


(3.119) 


V Z n< * Zr 7 avgy(JJ) tj” 2 / ; ; , 1\ u” 2 d ,* I\ 

- -——Hz (ij + J-Hz (tj ~ ~) 

J 1 Vq Zq ■ A/ • <TgygyV ? yj L 

2 V grid Z r,avgy(j?f) 

| 1 Vp Zq-A/- Gavgxi}?]) 

2 V W Sr,avgx(i,j ) pn- 1 /.- A 

^—77- 7~~K' h y vJ) 

J I 1 v 0 Z 0 • A/ • Gqvgxjlj) 

2 ^ £ raV gX(i,j) 

ry V 0 1 


(3.120) 


°v^e (V) + 1 ^ 

1 Vq ^0 * A/ ' CJ avgx\hj) L 


1 + ^—^-“-^ 

2, Vgrid £r,avgx(j,j) 


a H”~ l (i,j) - To 


vo 1 

V 8 nd \Xr,avg(j,j ) 


(3.121) 


(/• j - f) - o j+ io+o+1 j) - o - i j) 


The TE electric and magnetic field grid update equations for non-conductive 


media (o=0)are 


ExQJ) » E " 1 O'j) + z ° v g °nd er d(i'jj ' Hz 2 ( ij + ^~ Hz (3.122) 

E n y (jJ) "Ef'QJ) +Z o^ £r ^ (/J) • [^" 4 (/- ij)-^""(/ + i,/)] (3.123) 


(3.123) 
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wj)=«r'oj)+ 


V gnd \X.r,avg(},j) 


(3.124) 


Ef\iJ+ 1) - E7\i,j - i) +Ef\i ~\J)-E7\i + \,j) 


F inally the free-space TE Z grid equations are 


EWJ) *E» x - l (iJ)+Zo^ • My + i)-ft 2 (V-1> 

K(U)* £r'(v)+Zo^ • [ + \j> 


(3.125) 


(3.126) 


(3.127) 


Hr\U) + Yo^ • E7\iJ + \)-E7\iJ-\)+Ef\i-\j)-E7\i + \j) 


Dual TM Z grid equations for conducting media are 


(312*) 

+ + (3-129) 


j 1 Vp Zq • <7 avg (/,_/) * A/ 

2 Zravgji^j) 

1 1 VQ Zq ■QavgjiJ) ■ A/ 

2 V gnd £ r,avg(j,j) 


Er\iJ)+ 


(3.130) 


v 0 1 


j _j_ 1 vp *- -* 

2 ^grid £, r avg(j' 1 j} 

The TM Z electric field grid equation simplifies for non-conductive media 


E"(i,j) * E n z 1 (i,j)+Zo 


vp 1 

Vgrid &r,avg(j,j) 


(3.131) 


Hx{i,j- \)~Hx 2 (iJ + \)+Hy 2 {i + \,j)-Hy 


Finally, the free-space TMz grid equations are 

ffioj)= ht 1 o j. )+■ UrW -1) - sT* o'j+1) 


(3.132) 
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H n y (i,j) = + Yo ^:• E z \i + \,j)-E z ) 


(3.133) 


E" z (iJ)*Er\i,j)+Zo^- d - . 
tf” 4 (/J -1) - tf” 4 (4/ + i) +/£*(' + ~ \J) 


(3.134) 


Note that TE Z electric and magnetic field grid equations in 2-D have the following 
general form 

E n r = Cex\E° u + Cex 2 VHf d (3.135) 

E™ = C Ey xEf + C E yiVHf (3.136) 

H n z ew = CmHf + Ctm^Ef + C m VEf (3.137) 

where C's are constants (real numbers) that depend on the media properties and the 

velocity ratio v</v - d , and "del" operator represents the spatial derivative (gradient). 


S imilar ly the TM Z electric and magnetic field equations can be written as 

Hr=c Hxl m ,d +c H *vE° id 

Hr = C H yxHf + CnyiVEf 
E n r = C E iE° z ld + C Ex 2 VH 0 x ld + C E yiVHy ld 


(3.138) 

(3.139) 

(3.140) 


2. 2-D Grid Termination 


The grid equations derived previously are valid for all the nodes except for the 
nodes on the grid edges. The reason is that a grid edge node for a field component in a 
transverse-to-z plane has only one neighbor node instead of four like a node that is not on 
the grid edge, as shown below. 
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boundary 


Figure 3.12. Edge Nodes of a 2-D Grid. 

An equation for the circulation of the z-directed field component can not be 
written for an edge node because a z-directed field component on the edge shown in the 
figure above can not be updated using the grid equations discussed so far. Extending the 
ideas of transparent grid termination (TGT) and the discrete boundary impulse response 
(DBIR) applied for 1-D problems, the update equations for edge nodes can be devised 
based upon convolutions of the fields one layer inward and the pre-calculated impulse 
responses from the inward layer to the grid edge nodes. An edge field will be expressed as 
a superposition of responses, as illustrated in the figure below for a top edge grid node. 
(The same applies to other edge nodes as weDL) 
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Figure 3.13. Edge Field as a Superposition of Responses. 

The superposition can be also "visualized" using the concept of a multiport. The 
nodes on grid edges can be considered as output ports of a linear multi-port. Assuming, 
for simplicity, that the grid is square with N edge field nodes on each side, the total 
number of output ports will be 4(N-1). The nodes on the layer just inside the grid 
boundary will be considered as input ports of the linear multiport. There will be a total of 
4(N-3) input ports. The equivalent multiport is shown below. 



Figure 3.14. Modeling 2-D Grid Boundary as a Multiport. 
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The input and output fields can be either electric or magnetic fields. Furthermore, 
the input and output fields need not be of the same "kind", that is the input fields could be, 
in principle, E-fields and the output fields could be H-fields, or vice versa. However, we 
will select the input and output fields to be of the same "kind", that is either all E or all H 
fields, such that the impulse responses h^t) will be "dimensionless". The input and 
output fields will be H-fields for TE Z and E-fields for TM Z fields. In this manner, the same 
impulse responses h mjj (t) can be used for both TE Z and TM Z grids. The equations will be 
derived for TE Z fields, that is with electric fields at the grid edges. Dual equations for the 
TM Z fields can be obtained simply by replacing E's with Hs in the TE Z equations. 

There are four sides of the square grid boundary. The top side has the electric 
field E z (0,y,t), the bottom side has the electric field E z (L,y,t), the left side has the electric 
field E z (x,0,t), and the right side has the electric field E z (x,L,t). The layer just inside will 

have the corresponding fields as E z (Al,y,t), E Z (L-Al,y,t), E z (x,Al,t), and E z (x,L-Al,t). The 

limit s on the values of x and y for the edges are 0 and L (0 < x,y < L), while the limits for 
the layer just inside are A1 and L-Al (A/ < x,y <L- A/). Using the grid coordintes (ij) the 

fields on the grid boundary will be, starting at the upper left comer (coordinates 0,0) and 

using a clockwise reference direction£'"(0,y), E”(i, N— 1 ),E"(N- l,j),E"(i,0). These can 

be grouped into a "vector" of electric fields on the grid boundary 

E boundaty (t) = [E z (0J,t) E z (i,N-l,t) E Z (N- l,j,t) E z (i, 0,/)f (3.141) 

where T indicates transposition (the vector is a column vector, but it has been written as a 

transpose of a row vector to save space). In an analogous manner an electric field vector 

can be formed for the layer just inside the grid boundary 

E jusUnside (t) = [E z (lJ,t) E z (i,N-2,t) E Z (N-2 j, t) E z (i,l,t)] T (3.142) 
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The advantage of the vector arrangement is that we can now use a single index 
(say index m) to identify a field on the grid boundary and another single index (say index 
n) to identify an electric field on the layer just inside the grid boundary. The use of single 
indices m and n allows for a direct and compact refefrence to the multiport representation 
of the grid termination problem The electric field on the grid boundary E^^t) can be 
expressed as a convolution of the electric fields at the layer just inside the grid boundary 

E just ^^(t) and a matrix of impulse responses H(t). 

Eboundaryif) ^f(0 * ^jWiT Jnsideif) (3-143) 

where the 4(N-1) by 4(N-3) matrix of impulse responses can be written as 

h\z(i) 1,3(0. h 1, 4 (A'- 3)(0 

h 2 ,l(t) . ^2,4<W-3)(0 

H(t) = . (3 144) 


[ ^4(V-l),l(0 h 4( N-l)a(0 .),4(V-3) (0 J 

The electric field at a botmdary node (identified by a particular value of the index 
m) can be written as the product of the m-th row of the impulse response matrix H(t) and 

the col umn vector of the electric fields just inside the boundary (indentified by index n) 

Eboundaryim^t) — (3.145) 

S ^" 3) /w(0 * EjustJnside(n, 0 = 1 ”^ V_3) {' Ejustjnside(p~> t ) ' h m ^(t - x)dx 
where the summation is over all the nodes on the layer just inside the boundary. The 

discretized forms of the above equations, involving samples of electric fields at t=kAt, are 
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the sums of discrete convolutions. Since a discrete convolution is itself a sum, the result 
for the boundary nodes will involve double summations (the expression below uses 
superscripts for time samples and convolution summation). 

*L**<»>) = (3.146) 

The convolution equations express the fields at the edges as weighted sums (the 

weighting coefficents are the values of the discrete boundary impulse responses h^) of the 
time histories of the fields "just inside" the grid. In that respect they are equations of the 
same type as the equations for the non-edge grid nodes (equations ?? and ??), except that 
they involve, in general, s umma tions with many more terms. However, the impulse 
responses h irM1 (t) converge to zero relatively fast which effectively reduces the number of 
significant terms in the convolution sums. The convolution sums can thus be truncated 
such that only a certain number of the most recent values of the electric fields just inside 
the boundary are needed (this reduces the number of electric fields that have to be stored 
and updated). The discrete boundary impulse responses h^kAt) need to be determined 
(and saved) only once, just like in 1-D, for a selected grid velocity v^^Al/At. The 
impulse responses can be stored in a rectangular matrix form. A particularly suitable 
matrix form would have 4(N-1) rows, one row per boundary node, and 4(N-3)N h 
columns, where N h denotes the "truncated" length of the discrete boundary impulse 
responses. This form is suitable because the boundary electric fields can then be 
calculated as a product of such a matrix and a column vector constructed of 4(N-3) 
subvectors of length N h . The subvectors represent the N h most recent values of the 
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electric field at nodes just inside the boundary. (The choice of N h depends on the desired 
amplitude "resolution", with the values over 20 giving very good results for v grid /v 0 = Jd 
, D-2 for two dimensional problems). These subvectors are updated at each time step by 

shifting all of their values down (that is "futher into the past") by one, discarding the last 

value, and updating the first subvector element with the most recent field value calculated 

via standard grid update equations. The double summations are thus efficiently executed 

as matrix-vector multiplications. 

The discretized impulse responses br^ are determined in a manner analogous to 
the 1-D case, using the discrete equivalent of the delta function which we will denote as 
5 k . It is convinent to use the equivalent multiport to explain the procedure. In order to 
determine the discrete boundary impulse response for a particular input port n, all input 
ports except the n-th input port are set to zero and delta impulse is applied to the n-th 
input port. Recording all the outputs from t=0 to t=N b At provides us with 4N discrete 
boundary impulses of duration N h . 

= (3.147) 

= 0 for p=\...4(N-3) and p*n (3.148) 

The process is repeated 4(N-3) times (since there are 4(N-3) input ports) and the 

results are stored in a 4(N-1) by 4(N-3)N h matrix. This matrix is then used to implement 
the transparent grid termination via the matrix multiplication by the 4(N-3)N b column 
vector of the time-histories of the fields just inside the boundary. Depending on the values 
of N h and N, there will be input nodes that are sufficiently for from an output node such 
that the time it takes for an input impulse to propagate to the ouput port exceeds N h . In 
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such a case, the contribution of such an input node to the output node is known to be zero 
in advance. This may be used to reduce the number of operations in a TGT 
implementation. The process of obt ainin g the 2-D discrete boundary impulse responses is 
carried out on a "large" grid whose size should be at least (N+N h by N+N h ), or 
equivalently, the distance from the TGT boundary to the boundary of the "large" grid 
should be at least N h Al/2. The termination of this "larger" grid (in typical applications 
N»N h and the "larger" grid would be only incrementaly larger) is immaterial, since the 
reflections off its boundary do not arrive at the output nodes (where they would have 
corrupted the impulse responses) before the time-stepping has been terminated. When 
calculating the discrete boundary impulse responses (DBJR), the nodes interior to the 
TGT boundary need not be updated since all the nodes on the layer just inside the TGT 
boundary (the input ports) are zero for t>0 (at t=0 only one node on the layer just inside 
the TGT boundary has the value of 1). The DBIR calculations thus require updates only 
for the nodes between the TGT boundary and the large grid boundary which can result in 
significant computational time savings if N»N h . The update equations used for the 
region between the "inner" (TGT) and the outer boundaries are the free space equations, 
since this region is assumed to be free space. (The medium between the two boundaries 
can be any homogenous medium extending to infinity, but the most practical case is the 
free space.) 
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field updates uuly lufiie shaded raffia: 



Figure 3.15. 2-D Grid for Calculation of TGT Boundary Impulse Responses. 

The shape of the discrete boundary impulse response h^t) determined as 
described in this section will depend on the relative positions of the observation (output, 
index m) and the source (input, index n) nodes. In general, the more the indices differ 

from each other the following will be noticed 

• the impulse response will start later; 

• the ma ximum value of the impulse response will be smaller; 

• the impulse response will be more "spread out" in time. 

To illustrate these, three sample impulse responses are shown. The observation 
point (output port) is the same for all three responses, and is positioned in the center of 
one side of the TGT boundary. The three discrete boundary impulse responses are the 
responses to the delta function applied to 





• a source (input node) immediateUy below the output node; 

• an input node 5 nodes to the side of the source node immediately below the 
output node; 

• an input node 10 nodes to the side of the source node immediately below the 
output node. 

The velocity ratio for the DBIR calculations was set to 1/J2 . The plotting 
program interpolates linearly between the sampling points kAt. The impulse response 

were truncated to a duration of N h =40. Because of the vast change in scale (the largest 

value of DBIR for a source node just below is 1/2 while the largest DBIR value for a 

source 5 nodes to the side is about 25 times smaller, and the largest DBIR value for a 

source node 10 nodes to the side is about 1,000 times smaller) the impulse responses are 

shown individually. The rapid decrease of the impulse response peak values with distance 

offers the possibility to implement TGT using "local" equations instead of using "global" 

convolutions. An approximate "local" TGT implementation would involve only a few 

input nodes nearest to an output node, while a "global" TGT implementation involves all 

input nodes (all the nodes on the layer just inside TGT). The loss of accuracy due to 

"localization" is not excessive, because of the rapid reduction of the maxima of h^t) as 

the difference between m and n increases. Finally, the DBIR shown are equally applicable 

to TE z and TM Z grids. 
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Figure 3.17b. DBIR For A Source Node 10 Nodes to the Side 

C. 2-D TGT RESULTS 

Li order to test the TGT we will use the TMj grid model shown in Figure 3.11. 
Therefore, the nodes on the grid edges are E-field nodes. Note that the results would he 
the same for TE Z fields, with H-fields node on the edges of the grid. The source is a 
unit-amplitude delta impulse of E z electric field applied at the center node. The impulse 
source at the center of the grid gives rise to a cylindrical wave that propagates radially 
away from the grid's center. The cylindrical wave is not perfect, as its is "distorted" by the 
grid and develops a "wake" as it propagates through the grid. The cylindrical wave is 
incident to the TGT boundary. If the TGT were perfect, the total power within the grid 
would be the same as the total power within the same region for an infinite grid. 
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However, since the TGT is not perfect (because of truncaton of impulse responses) there 

will be some increase in the total power that comes from the "reflection" off TGT. A 

measure of TGT performance in 2-D will be the relative increase in the total power, that is 

the ratio of the power increase for the TGT and the power for an infinite grid. This ratio 

is a function of time and, since it is normally very small it is best represented in dB 

^ , (Pin[-PTGT) 

Qtgt = 10 • log 10 ^—J 

We have denoted this ratio as Q TGr since it expresses the TGT "quality". We will 
select a square grid with 21 nodes on each side for our 2-D TGT test. The impulse 
responses will be truncated such that the TGT performance can be observed for various 
selected impulse response durations. Note that the power within the grid will be the same 
for the case of an infinite grid and the TGT until the outward-propagating cylindrical wave 
reaches the TGT. All results will be shown for the case 

V grid = 1/2" • Vo 

The grid power and the TGT quality factor Q TGr will be plotted as functions of 
time, with the coordinate origin (t=0) reffering to the time when the outgoing cylindrical 
wave first reaches the TGT. 

First, we truncate the impulse responses such that their duration is t h = lOAt. The 
power within the grid is plotted as well as the power difference and the relative power 
difference in dB (the Q TCT ). The plot of power vs time shows the decrease in power 
towards 0 as the cylindrical wave leaves the grid. If the grid and the TGT were ideal, the 
power will be constant until the wave reaches the grid nodes at the centers of each side 
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(the outgoing wave reaches these nodes first since they are closest to the grid center 
where the wave has originated). The power will then decrease towards zero as more and 
more of the wave leaves the grid. Finally, the wave would reach the four comers of the 
grid (these nodes are the farthest from the grid center and the wave gets to them last) and 
leave the grid alltogether, the power within the grid being zero from then on. However, 
since neither the grid nor the TGT are perfect, there will be some small "residual" energy 
within the grid at all times, converging to zero as time progesses. The relative power 
difference in db (the Q TCT ) curve shows that, for t< t^ the power "reflected" off the TGT is 
about 150 dB below the power that the same domain (within the TGT boundary) would 
have for an infinite grid. For t > t h the power "reflected" off the TGT is more that 20 dB 
below the power of an infinite grid. This means that the effects of the TGT are about two 
orders of magnitude smaller than those of the grid "imperfection" for t > t h and about 15 
orders of magnitude smaller for t < t^. The plots for t h =20, t h =30, and t h =40 confirm the 
above as well Note that increasing the duration of the impulse response also reduces the 
"reflections" off TGT for t>t h . Also, note that for a grid of N nodes on the side the 
duration of impulse responses greater than 2N allows for any node on the layer just inside 
TGT to contribute at least one value to any of the nodes on the TGT (because there are 
2N nodes from one comer of the TGT to the opposite comer). 
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Figure 3.18. Residual Power for t h =3At. 
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Figure 3.19. Residual Power Difference Between TGT and Infinte Boundary for t h -3At. 
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Figure 3.20. Residual Power for ^=10At. 
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Figure 3.22. Residual Power for —20At. 
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Figure 3.28. Maximum Power Difference Between TGT and Inflate Boundary as a 
function of Impulse Response Duration (t h ). 
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Residual Power for TGT Boundary for vr=1/[2*sqrt(2)] 
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Figure 3.30. Residual Power for t h =40At, using Vr = 
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Figure 3.31. Residual Power Difference Between TGT and Infinte Boundary for t^OAt, 

using Vr = . 
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IV. ANALYSIS OF TGT FOR 3-D FD-TD 


A. FD-TD FORMULATION IN 3-D 

1. FD-TD Equations in 3-D 

The incident and the scattered electromagnetic fields and the media parameters in 
3-D problems vary with three spatial coordinates(x,y,z) , so the electric and magnetic 
fields will be functions of the spatial coordinates x, y, z, and time t. A 3-D electric field 
vector E can be represented in terms of its othogonal components as shown below 


A z 


Ex 


.-.» y 


Figure 4.1. Electric Field Vector Components. 

~E(x,y,z, t) = E x (x,y,z, t)f +E y (x,y,z, t)~y + E z (x,y,z, t) z (4.1) 

The direction of propagation is indicated by the direction of the velocity of 

propagation vector ~v. Similarly, for the magnetic field vector H we can write 

H(x,y,z,t) = H x (x,y,z,t)x +H y (x,y,z,tjy +H z (x,y,z,t) z (4.2) 

The electric and magnetic fields are solutions of Maxwell's equations 
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(4.3) 


f Ct E(x,y,z,t)'!u E = -|{Jf s , • **} = -|{JL, vix,y,z)H(x,y,z,i)-lZ 


» -» 


| Cftr 0 • din = 


(4.4) 


a. 

dt 


Jls w <x,y,z)E(?,y,z,t) • dS • + o(x,y,z)E(x,y,z,t) • 


The second curl equation does not have the source current term, since it has been 
assumed that there were no source currents in the domain of interest. The contours of 
integration for the electric and the magnetic field circulations are, in general, different and 
will thus labeled C E for an electric field contour and C H for a magnetic field contour. 
Similarly, the surfaces associated with the contours will be labeled S E for the magnetic flux 
and S„ for the electric flux. In order to determine the discretized forms of Maxwell's curl 
equations in the integral form for the 3-D, we will use 3-D Yee electric and magnetic cells 
[Ref I]. One such cell is shown below. 
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The electric field contours C E and the magnetic field contours C H are square (A1 by 
Al) contours. The linked C E and C H contours are all orthogonal to each other, evoking the 
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idea of electric and magnetic field linkage in all directions. Based on Yee's cells, the 
electric and magnetic contours and their associated surfaces will be used to discretize 
Maxwell's curl equations. 

2. Derivation of Magnetic Field Update Equations 

We need to derive the update equations for each of the three components of 
magnetic field vector H. A contour C E for deriving the update equation for H^y^t) is 
shown below. This contour can be thought of as the front face of the Yee cell (cube) from 
the previous figure. Note that the coordinates of the four electric field components are 
expressed in terms of the coordinates of the contour's center. 











The local coordinates will be used in the evaluation of the line and surface integrals 

that constitute the integral forms of Maxwell's equations. The circulation of the electric 

field around the electric field contour C E , assuming counter-clockwise reference direction 

such that the normal to the surface S E is in the direction of the magnetic field H x (x,y,z, t), 
can be evaluated approximately, assuming that the electric fields are constant over their 

respective countour's sides of length A1 

jy E y (x',y',z',t)x +E z (x',y',z',t)~y ] • dl E « (4.5) 

E y (x,y, z-y,0~ E y (x,y,z+ ^,t) +E z {x,y + y, z, t) - E z (x,y - y, z, t) ■ M 
The magnetic flux through the surface S E (the right-hand side of the first curl equation) 
can be calculated using the local coordinates 

jj Si n(x / ,/,z / )// z (x / 5 /,z / , t)~z • ds E = (4.6) 

J J J J \i(x,y + q/,z + q)H x (x,y + \\t,z+q,t)x * xdydq 

2 2 

There are infinitely many ways to model the variation of the magnetic flux density 

(with respect to local coordinates vp and Q within the contour. We first need to postulate 

a certain variation of the magnetic field with vp and C, such that the integral can be 

evaluated over S E . The simplest model assumes that the contour width A1 is small enough 

such that the magnetic field H(x,y+vp,z+£,t) within the contour may be assumed constant 

and equal to the value at the contour’s center H(x,y,z,t). 

H x (pc,y + vp, z+q, t) « H x (x,y, z, t) (4.7) 

This is equivalent to using a piece-wise constant (or 2-D "pulse" expansion) 

approximation of the actual magnetic field variation with respect to y and z. Note that the 
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above assumption allows that, although the magnetic field is constant within a contour, it 
can change from one contour to an other. This yields an approximate expression for the 
for the magnetic flux 

ff v(x> ,y',z')H x {x',y',z',t)~z • ds s *H x (x,y,z, 0 • fj \ + l \x(x,y+y,z+q)d\ydq (4.8) 

JJ Ss 2 2 

The integral of the permeability p(xy+v|/,z+Q can be re-written in the following manner 
f! f| v(x,y + vp, z + q)cTqd\y = (A/) 2 • 7777 fl fl z + ( 4 - 9 ) 

j _ j —— IA /1 •> •> 


,y + W, Z + q)d\ydq (4.9) 


The term in the brackets is recognized as the permeability averaged in the y 
and z-directions within the contour C E (which is in a x = const plane) 

tW(x,y, z) = —77 J + | \ + _l \i(x,y + 1 |/, z + q)d\ydq (4.10) 

(A/) 2 2 

The approximate expression for the magnetic flux through S E can now be written using the 
average permeability as 

JJ p(x',/, z')H x (x',y f ,z',t)~x *ds E ~ (Af) 2 [Lav g x(x,y,z)H x (x,y,z,t) (4.11) 
Ag ain this approximate expression resulted from the piece-wise constant 
approximation of the magnetic field with respect to y and z-coordinates. The main 
advantage of the piece-wise constant expansion employed above is its simplicity. The first 

curl equation of Maxwell can now be replaced by an approximate equation 

E y (x,y, z - y,0- E y (x,y, z + yt) + E z (x,y + y, z, t) - E z (x,y - y, z, t) * (4.12) 
-y{ [\i.avgx(x,y, z)H x (x,y, z, t)]A/> 

which can be simplified (because of media stationarity) to 

E y (x,y, z — y, f) — E y (x,y, z + f,t)+E z (x,y + z, t) - E z (x,y - z, t) « (4.13) 

-^■{H x {x,y,z,t)} ■ Al■ p mgx (x,y,z) 
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This equation can also be re-written as 




(4.14) 


77 - - -- \e,(x, y,z - /) - E y (x,y,z + ^-,t)+ E z (x,y + y ,z, t) - E z (x,y - y, z, I) 

The time derivative operator in the above equation is typically replaced by the 
finite difference approximation [Ref 5], However, just like in 1-D and 2-D, we present an 
alternate approach, such that the approximation of the field time variation is shown to be 
analogous to the approximations already introduced for the field spatial variation. The 
above equation is integrated with respect to time to get an approximate equation for the 


magnetic field at the present time t 


H x (x,y, z, 0~ A/ . z ) x 


(4.15) 


[f 0 E y (x,y,Z - f , x)dx - {' E y (x,y, z + f , x)dx + £ E s (x,y + f,z, x)dx- {' E : (x,y - f , z, z)dx\ 

A similar integral equation can be written for the magnetic field at the same spatial 


location but at an earlier time t-At 


H x (x,y, z,t- At) a —----x 

Al • iXavgx(.X,y,z) 


(4.16) 


J' ^ Ey(x,y,z- y,x>fx-J' * Ey(x,y,z + f ,x)dx + J' * E : (x,y + ^,z,x>fx- J' * E z (x,y- f ,z,x>/x] 


Subtracting the two equations we get 


H x (x,y,z,t)aH x (x,y,t-At) - —----x 

Al ■ Havgx(x,y,z) 


(4.17) 


Ey(x,y, z-f, x)dx - Ey(x,y, z + f, x)dx + E z (x,y +f,z,x)dx- E : (x,y - z , x>/x] 


The integrals on the right-hand side can not be evaluated exactly, because the 
exact temporal variation of the electric fields within the At interval prior to t is generally 
not known (just like the spatial integral for the magnetic flux could not have been 
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evaluated exactly because the exact variation of the magnetic field within the contour was 
not known). However, the integrals can be evaluated approximately by assuming a certain 
variation of the electric field with the temporal variable t. The simplest approach, 
consistent with the assumptions made for the field spatial variation, would be to assume 
that the interval At is small enough such that the electric field can be assumed 
approximately constant within At and equal to the value at the center of the tune interval 


(t-At, t) 


Ey(x,y,Z~ y,T)s 

*E y (x,y,z-j,t-j) 

for 

t-At <x<t 

(4.18) 

Ey(x,y,z+ y,T)s 

aE y (x,y,z + j,t-f) 

for 

t-At<x< t 

(4.19) 

E z (x,y + f,z,x)< 

*E z (x,y+j,z,t-j) 

for 

t- At <x<t 

(4.20) 

E z (x,y-f,z,x)t 

*E z (x,y-f,z,t-f) 

for 

t-At<x<t 

(4.21) 


The above represents a piece-wise constant approximation of the electric field 
variation with respect to the temporal variable t. The approximate expression for the 


magnetic field at location (x,y,z) and at the time t now becomes 

H,(x,y,I) ««,<*,*^' (w j » < 4 - 22 > 

[E y (x,y,z-f,t-f)-E y (x,y,z+f,t-%)+E z (x,y+f,z,t-f)-E z (x,y-f,z,t-f)] 
The above equation can be written in a somewhat different form, using the identity 

1 1 1 _ J&0 l_l_1_1_1 __L_y_L (aox 


i _ j_j_ i i = i i_ l = v ii = voTo— 

V HO Vr ys °Z 0 ^ ^ 


(4.23) 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Z 0 * 377H), and p r denotes relative permittivity. Introducing the 
grid "propagation" velocity 

< 4 ' 24 > 
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we can write the approximate "update" equation for the magnetic field as 

A/) - Fo(^j) * ( 4 - 2: 

[£,(*> y,z -f,t-f)-E,(x,y,z+f,t-f)+E z (x,y+f,z,l-j)- E y (x,y - f ,z, 1 - f ] 


(4.25) 


The equation simplifies for the case of non-magnetic media (p=l) 


H x (x,y,z, t) « tf,(x,y, 2 ,7-A/) - To-J x 


(4.26) 


0,(*,.y, z - f, t - f) - E y (x,y, z +f. t - f) *E z (x,y +f ,z, t - f) - E,(x,y - f, z, t - f )] 

The update equations for the components Hy(x,y,z,t) and K, (x,y,z,t) can be 
derived in the same mann er. Since the only differences would be due to different contour 
and surface orientations, the update equations for the components Hy(x,y,z,t) and K, 
(x,y,z,t) can be obtained from the update equation for H x (x,y,z) by cyclic permitation of 
indices. The contour for updating Hy(x,y,z,t) is shown below 


A z 


Ez(x+dx/2,y,z,t) 



Ez(x-dL / 2 > y,z,t) 


Ex(x,y,z-dl/2,t) 


Figure 4.4. A Contour for Updating 


Following the same procedure shown for updating component, or simply using 
the cyclic permutation of indices we get the update equation for Hy(x,y,z,t) 
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H,(x,y,z,t)*H y (x,y,z,l-At) - r„Q^) * < 4 - 27 ) 

[E z (x - f,y, z, t - f )-E z (x + f ,y,z, t - f ) +E x (x,y,z + E x (x,y,z- f , t - f )] 

Finally, the contour for updating H z (x,y,z,t) is shown below 


Ex(x,y-dl/2,z,t) 


Ey(x-dl/2,y,z,t) 


Hz(x,y,z,t) 


Ex(x,y+dl/2,z,t) 


Ey(x+dl/2,y,z,t) 


Figure 4.5. A Contour for Updating H,. 

The update equation for is 

H,(x,y,z,i)*H,(x,y,x,t-A,) - Yo(^-) * (4.28) 

[E,(x,y-f,z,t-^-E x (x,y + f,z,t-f)+E y (x+f,y,z,t-f)-E,(x-f,y,x,t-f)] 
In most cases the media will be non-magnetic, in which case all relative 

permittivity averages will be equal to 1, and the update equations will simplify accordingly. 
3. Derivation of Electric Field Update Equations 

We will start with the Maxwell's curl equation for the circulation of the magnetic 
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j c ^H z (x,y,z,t)~z * dl H = j f jjj^ e (x,y,z)^E x (x,y,z,t)x +E y (x,y,z,t)'y]* ds H +(4.29) 

a(?C ’ y ’*)x +E y (x,y,t)~y ] • ds H 
The second term on the right-hand-side (the conduction current) is obviously 
going to cause the update equations for the electric field to he more comphiicated than the 
magnetic field update equations. The discrete equivalent of the magnetic field circulation 
needs to be determined for three sample contours a contour in a plane parallel to the 
yz-coordinate plane, a contour in a plane parallel to xz-coordinate plane, and a contour in 
a plane parallel to the xy-plane. The first contour provides an update equation for the 
x-directed component of the electric field, the second contour provides an update equation 
for the y-directed component of the electric field, and the third contour provides an update 
for the z-directed component of the electric field. Again, it is sufficient to derive the 
update equation for one electric field component only (say the x-component) since the 
other equations can be obtained by a simple cyclic permutation of indices. 

Z A 

i 

Hz(x,y-dy/2,z,t) 

^ Hy(x,y,z-dz/2,t) 

Figure 4.6. A Contour for Updating E x . 
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A local coordinate system (\|/, <;) will be used, with the origin at the center of the contour. 
Any point (y^z 1 ) within or on the contour C H (or, equivalently, any point on the surface S^) 

can be specified by its local coordinates 

y'=y + \ j/ z‘ = z + q where - <q< +^- , and -^<v|/<+^ 

The circulation of the magnetic field around C H , assuming counter-clockwise 

reference direction such that the normal to the surface S H is in the direction of the electric 
field E x (x,y,z, t ), is approximately 

[H y (x'J,z', t)~y +H z {x , ,y',z', ifz ] • dl H « (4.30) 

[H y (x,y,z-f,t)-Hy(x,y,z + f,t)+H z (x,y + f,z,t)-H z (x,y-%,z,t)] -A/ 

The electric flux through the surface S H can only be evaluated approximately , since the 

exact way that the electric flux density 

D x (x,y, z, t) = s(x,y, z)E x (x,y,z, t) (4.31) 

varies with y and z is not known a priori. First we re-write the flux integrals 

J* eC*',/,*', t)E x (x'y,z', t) x *ds H = (4.32) 

JJ J _l s (x,y + V,z + q)E x (x,y + \\>,z+q,t) x • x d\\idq 
2 2 

JL„ q( x ' y i z ') e *( x/ y > z ' > t)t •ds H = ( 4 . 33 ) 

r+~ r+— —^ ^ 

\J, JJ a(x,y + \v,z + q)E x (x,y + \v,z+q,t)x • xdydq 
2 2 

Next, we assume that the contour side A1 is small enough such that the electric 
field E(x,y+v|/,z+£,t) within the contour may be assumed constant and equal to the value at 
the contour's center E(x,y,z,t). This assumption yields approximate expressions for the 
flux integrals 

JJ^ s(x,y + y,z + q)E x (x,y + y,z+q,t^x (4.34) 

• ds H « E x (x,y, z,/)-JJ J J z(x,y + \\),z + q)dydq 

2 2 
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jj^ o(x,y + \y,z +q)E x (x,y + vj/ 3 z+q, t)~x (4.35) 

•ds H ~ E y (x,y, J J, cs(x,y + q/, z+ <q)d\\/dq 

2 2 

or, using the permittivity and the conductivity averaged in the y and z-directions 
J} z(x,y + \y,z + q)E x (x,y + y,z + q, t)~? •ds H * (A/) 2 'E x {x,y,z,t)-z mgx (x,y,z) (4.36) 
]f a(x,y + q/,z + <;)E x (x,y + y,z+q,t)x *ds H a (A/) 2 ■E x (x,y,z,t)-c avgx (x,y,z) (4.37) 

The second curl equation of Maxwell can now be replaced by an approximate equation 

[H y (x,y,z-j)-H y (x,y,z + f) +H z (x,y + f,z,t)-H z (x,y-f,z,t)] -A/ *(4.38) 

j t {[Zavgx(x,y,z)E y (x,y,zJ)](kt) 2 } + (A/) 2 -a^(x,y,z) -E x (x,y,z,t) 

The approximate curl equation can be re-written as 

§- t {E,(x,y,z,f)} (4.39) 

—- - - -\h,(x,v, z - f, i) - Hy(x,y, Z+J,t) +H,(x,y +f ,z, t) - H,(x,y - f,z, ()] 

A/ * Zavgx[X,y y Z) 

The above equation is integrated with respect to time to get an approximate equation for 
the electric field at the present time t 


E x (x,y,z,t ): 


(4.40) 


n x \*,y,*, V ~ 77” ✓ „ , *■ v*-™/ 

A/ * ScrygxV^, J 7 ,-2^) 

[J' /^(x,y, z - f, - J' tf,,(x,y, z + f, x)dx +J' H z (x,y +f ,z, x)dx - £ H z (x,y - f, z, x)*] 

_ 0 ^( W ) p 

z*vgx(x,y,z) J 0 v ' ^ 

Similar equation can be written for the electric field at the same spatial location but 


at an earlier time t-At 


E x (x,y, z,t- A t) ~ — ---- x 

M-Savgx(x,y,Z) 


(4.41) 


[j' Af H y (x,y, z-f, x)dx - J' * H y (x,y, z + f, x)dx + J' * H z (x,y+ f, z, x) - {' ^ //.(*,>- - f, z, T>*t] 

Gavgx(x,y,z) c t-&t 

J o toD* 
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Subtracting the two equations we get 


E,(x,y,z, 0 * E,(x,y, z, t-Al)- » 


(4.42) 


[J^ H,U,y,z- f.xyh- f'_ 4| H,(x,y, z+j,z)dz+^H ! (x,y+j,z,z)>h-jH 1 (x,y- f,z, x>*] 

- a -z- ( *-y- z ) I' E ,(x,y,z,x)dx 
£avgx(x,y,z) !'-* xK ^ 

Assuming that the interval At is small enough such that the magnetic fields can be 
assumed approximately constant within At and equal to the value at the center of the time 


interval At 


H y (x,y,z+f,x)*H.(x,y,z + for l-At<x<t 

Hy(x,y,z-f,x)«H y (x,y,z-f,t-f) for t-At<x<t 
H 2 (x,y + f,z,x)*Hfx,y-¥,z,t-f) for t-At<z<t 
H-Jx,y-f,z,x) » ff,(x,y-f,z,t- f) for t-At<z<t 


we get 


E x (x,y, z, t) « E x (x,y, z, t - At) + 


A/ • e^Cx^z) 


(4.43) 

(4.44) 

(4.45) 

(4.46) 


(4.47) 


[H y (x,y, z- f, t- f) - H y (x,y,z + +H z (x,y + f,z,t-f)~ H z (x,y -f,z,t- f )] 

_ o avgx (x,y,z) rr Ex(xy z/c)dT 

£avg X (x,y,z) t to 

The last integral term can be expressed in different form as shown below 


* - Stti -to ■ L- A ' • «- 48 > 

The above equation can be written in more compact manner by using the identity 

1 1 1 1 1_L fEi (4.49) 


1 — - V _ X I _ 1 / rv 1 _ 7 1 

S - s 0 s r -^-ye^s r - yp^-Vso S r 0 °s r 


where the free space velocity of propagation is denoted v 0 , the intrinsic impedance of free 
space is denoted Z 0 (Z 0 « 377Q) and e r denotes relative permittivity. Introducing the grid 
"propagation" velocity 


v 

Vgrid - A{ 


(4.50) 
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we can write the approximate equation for the electric field "updates" as 

E x (x,y,z,x)*E,(x,y,z,x) x (4.51) 

[>,(*, y,z - f , t - f ) ■ - HJx,y,z + f , t - f) +H,(x,y +f ,2,/■- 4 () - H,(x,y - f .2, t - f)] 


^ f>v gride ravgI (x y ') ^ 


The average of the electric field between t-At and t, denoted E x _ av , g(At) can not be 
calculated exactly because the exact temporal variation of the electric field within the At 
interval is not known a priori. However, assuming a linear variation within At interval, the 
temporal average of the field is the average of the field values at t-At and t 

(4.52) 

Substituting the time-average E V JVg(At) (x,y,t) into the update equation we get 

E x (x,y,z, t) * E x (pc,y,z, t - At) +Z 0 ^-- - j ——r x (4.53) 

v gndS ravgx (X,y,Z) 

\_H y (x,y,z -f,t- f) - H y (x,y, x + fj-f) +H z (x,y +f ,z,t- At) - H z (x,y -f,z,t- f)] 
^ v 0 1 a / _ . v x E x (x,y,z,t-At)+E x (x,y,z,t) 

' A ' <w(w) -2- 

which give the update equation for the x-component of the electric field as 

l 1 vq Z 0 ■Al-c avgx (x,y,z) 

E x (x,y , z, t) *- 2 Snd £ A 7f^; } v E x (x,y, z, t-At)+ (4.54) 

j [ 1 Vq • A/ • <?gvgj;(x,y, z) 

2 Vgnd s raV g X (x,y,z) 

Z — 1 

_ V X nd £-r,avgx(x,y, z) _ 

1 | 1 VQ Z 0 • AZ-a^fo^z) 

2 V £r^x(x,y,z) 

[/^(x,y, Z -f,t-f)-//>,y,x + f,t-f)+i/ t (x,y + f,z,t-Af)-i/ r (x,y-f,z,t-f)] 

The equation simplifies greatly for the case of non-conductive media (a= 0 ) 
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E x (x,y, z, t )« E x (x,y, z, t-At) + Z 0 


(4.55) 


^^V gndZravgx ( x ^ 2 ) ~ V""> 

[H y (x,y, z-f,t-f)- H y (x,y, z + f,t-f) +H z (x,y + j,z,t- At) - H z (x,y -f,z,t~j)] 
In the case of free-space (e =1), the equation simplifies further 


Vn 

E x (x,y,z, t) * E x (x,y,z, t- At) +Z 0 ^—- x 


(4.56) 


[H y (x,y, z-f,t-f)~ H y (x,y, x + f,t-f)+ H z (x,y+f,z,t- At) - H z (x,y -f,z,t- f)] 
The update equations for the other two electric field components E y (x,y,z,t) and E z 

(x,y,z,t) can be derived in the same manner. Since the only differences are in the contour 
and surface orientation, the update equations for E y (x,y,z,t) and E z (x,y,z,t) can be 
determined by a simple cyclic permutation of indices, as shown on the next page. The 
contour for updating E y (x,y,z,t) is shown below. 


Hz(x+dx/2,y,z,t) 



Hz(x-dI/2,y A t) 


Hx(x,y,z-dl/2,t) 


Figure 4.7. A Contour for Updating Ey. 

Following the same procedure as shown for the E x (x,y,z,t) update equation, or by 
a simple cyclic permutation of indices we obtain 
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E y (x,y,z,t) 


j 1 vp .Zq • A/ • Cqygy(X,y,Z) 

s (x,y) v _ Af) + (45?) 

2V^ £r,av2y(x,J,z) 


7 ' V _A_ 

° v gride ravgy (x,y,z) 

J ^ 1 Vq Zq • A/ • CTavgy(x,y,z) 

2 Vgnd z r>avgy (x,y,z) 

[H,(x- j,y,z, /-AO-#z(* + f,y,z,/-f) + H x (x,y,z + f, t- f)-tf,(x,.y,z-f,/-f)] 
The equation simplifies greatly for the case of non-conductive media (o=0) 


Vn 1 

^(r.y.2,0 « E y (x,y,z,t-M) +Z°^ w(] ^ z) * 


(4.58) 


[#z(x - f,y, z, t- At) - H z (x + J,y, z,t- f)+ H x (x,y,z +f, / - f) - tf*(x, j, z - f, / - f)] 
In the case of free-space (s r =l), the equation simplifies further 


Ey(x,ya, I) ~ E y (x,y,z, I - At) +Zo—;* 


(4.59) 


v — VgnJ v -> 

[H z (x - f ,y, z, t- A/) - tf z (x + f, y, z,t-f)+ H x (x,y,z +f, / - f) - H x (x,y, z-f,t- f)] 
The contour for updating E z (x,y,z,t) is shown below. 



Figure 4.8. A Contour for Updating E r 







Following the same procedure as shown for the E x (x,y,z,t) update equation, or by 
a simple cyclic permutation of indices we obtain 

j 1 Vo Zp * A/ • ggvgz(x,y,z) 

*#**-*>+ («■«> 

+ 2'W Sr^v«z(x,y,z) 

z v ° 1 

0 ^srids r<avgl (x,y,z) _ 
j ^ 1 Vo Zo • A/ - CTavg.(x,y, z) 

2 Vgrid S ravg! (X,y,Z) 

[H x (x,y - f ,z,t-f)~ H x (x,y + f, z, t -- f) + //,(* + f,y, z, / •- f) •- #y(* - j,y,z, t - f)] 

The equation simp lifies greatly for the case of non-conductive media (cf=0) 

E z (x,y,z,l)*E z (x,y,z,l- A/) < 46 » 

[ H.(x,y -f,z,t-f)-H,(x,y+f,z,l-f)+H,,(x+f,y,z,t-%)-H r (x - f,y,z, I -- f)] 
In the case of free-space (s =1), the equation simplifies further 


E z (x,y,z, t) « E z (x,y,z,t-Af) +Z 0 ^x 


4.62) 


[l H,(x,y f ) - H,(x, y + f,z,t-f)+ H r (x+f,y,z,t- f) - - fo’,*. /- f )] 

B. TRANSPARENT GRID TERMINATION IN 3-D 


1. 3-D Grid 

The electric and magnetic field update equations have been derived using local 
coordinate systems, with origins at the centers of the magnetic and electric contours, 
respectively. These equations now need to be "converted" to a global coordinate system, 
tha t is to a 3-D grid of equi-distant sampling points. We assume that our domain is an L 
by L by L cube discretized using a uniform grid, that is that the discretization step A1=L/N 
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is not chang ing with position. The electric and magnetic field sampling points (nodes) are 
"interleaved", that is offset by Al/2 from each other in x, y, and z-directions. The electric 
field components are located along the sides of "electric" cells while the magnetic field 
components are located along the sides of "magnetic" cells. 






The electric cells are the cells of constant s and a, while the magnetic cells are the 
cells of constant p. The permitivity and conductivity can vary from one electric cell to 
another in arbitrary manner. The permeabilty can vary from one magnetic cell to another in 
an arbitrary manner as well. The electric cells are shifted by Al/2 in all three coordinate 
directions with respect to corresponding magnetic cells. The placement of field 
components along the cell's sides automatically provides for the continuity of tangential 
components of electric fiels across media interfaces, while the spatial shift of electric vs. 
magnetic cells provides for electric and magnetic field linkage (coupling). Domain 
boundaries can coincide with either the faces of electric cells or the faces of magnetic cells, 
but not both (because of the Al/2 spatial shift between the two). Although, in principle, 
either electric or magnetic cells can be alligned with the grid boundaris, we will select 
electric cell alligment with the domain boundaries. 
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Figure 4.10. Yee Electric Cell and Associated Field Components. 

A grid of (N by N by N) electric cells has (N)(N)(N) E x , E r and E z electric field 
nodes, (N-1XN-1XN-2) H*, Hy, and H z magnetic field nodes (Just the opposite would hold 
for alligning magnetic cells with domain boundaries). 3-D electric and magnetic field 
"update" equations can now be "converted" to electric and magnetic field grid "update" 
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equations by replacing the variables x,y,z and t with the grid coordinates of the electric 
field spatial and temporal sampling points 

x -»/A/, / = 0 , 1...N y -» jAl, j = 0 , I...N z -» kAl, k = 0 , 1..N t -> «At, n = 0 ,1...#, 
The notation can be simplified by omitting the common A1 and At terms, and using 

a superscript for the index of the temporal sampling point. The E field grid equations for 

conductive media are obtained from equations (4.54,4.57,4.60) 

i *, 7 a . A / . rr (i i 1r\ 


E"(i,j,k) ■ 


1 _ 1 Vp Zo ‘ ^ * *3gvgxijiji k) 

2 *W Zr^vgxjiJ, k) rn-lfj • , 

n lV 0 Zo • A/ • GqvgxijJ, k) £ ' X } 

2 V grid g rjavgx (j,j, k) 


ry Vo 1 

^0 .- T. 


V grid £ r aV g X (i,j, k) 

l | 1 Vo Zo • A/ • g avgxjjj, k) 

2 V grid Sr4r»gx(x,y,z) 

h 7 1 (tj, k-\)~ H n ; h QJ, k+\)+H7 h (tj + \,k)~ e7 1 (/J -\,k 

1 _ j Vp ^0 • A/ ' <5avgy(i,j, k) 

E"y(!J,k) * A 777i iki E ^ iJ ' k) + 

] I 1 v 0 ^0 ‘A/ 'GqygyVjdi k) 

2 £ r f avgy (j J, fy 

■ rr Vo 1 

Zot;— -- r r":r 


(4.63) 


(4.64) 


_ £r,avgy(j>j, &) _ 

J ^ 1 v 0 • A/ * Gqygy(J,j? ff) 

2 V grid £r&vgy(ij, fy 

#r' o - \j, k) - h 7. i (i+\,j ,*>+#r ■ (/j, t+1) - <v, * - §> 
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(4.65) 


1 _ 1 vo Zo-U-^ ovgzjijfk) 

m,j,k )« f *" ^ (v.*)+ 

2 I 1 v 0 ^0 * ‘ Gqygz(jj> 

2 V grid Sr,avgz(iJ,k ) 


Sr l 0 ygz(/,y, At) 

J _|_ 1 v 0 ^0 ’ ' Gavgz(i,j->k') 

2 V grid S r,avgz(ij,k) 

h 7 2 (i,j~j,k)~Hx 2 (i,j+\,k)+HT 2 (i + ^,j,k)-Hy 2 (i-\,j,k) 
The equations simplify greatly for the case of non-conductive media (o=0) 


E n x {iJ, k) « E” 1 Q,j, k) + Z 0 


Vgrid 8 r,avgx(j,j, k) 


H7kiJ^-\)-H7kiJM\)+Hf 2 {iJ + \,k)-H7\i,j-\,k) 


Ey(iJ, k) « Ey 1 (i,j, k) + Zo 


V gnd Zrpvgy(},j, k) 


Hz 2 (/ - \,j, k)-H z 2 (/' + \,j, k) +H7 2 (i,j, k+\)~ h7 2 ( ij, k-\) 


(4.66) 


(4.67) 


E"Q,j, k) « E"~ x (ij, k) +Zo^-~ - ±—-x 

% nd Zr,avgz\l>J? k) 

Hx 2 (iJ-\,k)-H x 2 (ij + k) +H7 2 (i + \,j,k)-Hp(i- \,j,k) 
In the case of free-space (s=l), the equations simplify further 

EWJ,k)*ET\Uj,k) +Zo^x 
B7\iJ,k-\)-Hf'(iJ,k+\)+H?\i,j + \,k)-lfi i (i,j-\,k) 


h 7 ‘c - \,J 


E",(ij,k)*Er'(iJ,k) + Z 0 ^x 

j, k) - h7‘ (/ + U k) +h 7 1 (IJ, k + j) •- H? 1 ' ( ij, k-\) 


Et(iJ,k)»ET L (lJ,k ) +Zo^x 


h7' 1 (IJ - k)-HT‘(iJ + 5 , *) +h7‘ 0 + \j, k)- h7' (/■-± J,k) 


(4.68) 


(4.69) 


(4.70) 


(4.71) 


113 



Similarly, the H-field update grid equations are 


E n ; h (iJ , k-\)~ eT* (i,j, k + \) +Ez~\i,j +i A) -JST^U - *) 




yp 1 

V t"t> p.y(l,J,k) 




(4,72) 


(4.73) 


E7'(i,j-\,k)-E7'‘(i,j + \,k)+£" 1 (i + jJ,k)-Ey ! ( t-\j,k ) 


(4.74) 


2. 3-D Grid Termination 


The grid equations derived previously are valid for all the nodes except for the 
nodes on the domain boundary. The reason is that the field nodes on the boundary have at 
most five nearest neighbor nodes while the nodes within the domain boundaries have six 
nearest-neighbor nodes. 
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because the field component on the boundary shown in the figure above can not be 
updated using the grid equations discussed so far. Extending the ideas of transparent grid 
te rminat ion (TGT) the discrete boundary impulse response (DBIR) applied for 1-D and 
2-D problems. Similarity to the 2-D TGT concept, the update equations on the boundary 
can be devised based upon convolutions of the fields on the just inside boundary and the 
pre-calculated impulse responses from the just inside boundary to the boundary. The 
fields on the boundary will be expressed as a superposition of responses, as illustrated in 
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the figure below for top plane of the boundary. The same applies to other plane of 
boundary as well. 


the node on the boundary has 



Figure 4.12. Fields on the Boundary as a Superposition of Responses. 

The superposition can be also modelled using the concept of a multiport. The 
nodes on the grid boundary have three field components and those field components can 
be considered as output ports of a linear multi-port. Assuming that the grid is a cubic with 
(N)(NXN) nodes, the total number of output ports will be 3[2(NXN)+(N-2)4(N-1)]. The 
field components on the nodes just inside the grid boundary will be considered as input 
ports of the linear multiport. There will be a total of 3[2(N-2XN-2)+(N-4)4(N-3)] input 


116 





























ports. The equivalent multiport model in 2-D TGT can be used in the same manner. The 
input and output fields can be either electric or magntic fields. Furthermore, the input and 
output fields need not be of the same ''kind", that is the input fields could be, in principle. 
E-fields and the output fields could be H-fields, or vice versa. However, we will select the 
input and output fields to be of the same "kind", that is either all E or all H fields, such 
that the impulse responses h^t) will be dimensionless. The equations will be derived for 
E field model, that is with electric fields on the grid boundary. Dual equations for the 
magnetic fields model can be obtained simply by replacing E's with FTs. There are six 
planes for the cubic grid boundary. However, by 
makin g an arrangement for the input and output in 
TGT model as shown below. 


splitting the cubic into N pages and 
the N pages, we can simplify the 3-D 
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Figure 4.13. Page Arrangement for 3-D Cubic Cell. 

The boundary field components for the first page have the electric field E x (ij, l,t), 
E y (ij,l,t) and E z (Lj,l,t). The boundary field components for the second page through 
(N-l)th page have the electric field components like this. First, for the x components, the 
2nd Page has [E x (lj,2,t), E x (i,N,2,t), E x (Nj,2,t), E x (i,l,2,t)], nth (2<n<N-l) page has 
[E x (lj,n,t), E x (i,N,n,t), E x (Nj,n,t), E x (i,l,n,t)] and (N-l)th page has [E x (lj,N-l,t), 
E x (i,N,N-l,t), E x (Nj,N- l,t), E x (i,l,N-l,t)]. For simphcity, the y,z components have the 
same global coordinates used in the x components. Finally, the boundary field components 
on the Nth page have E x (ij,N,t), E y (ij,N,t) and E 2 (ij,N,t). The limits on the values of ij,k 
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on the boundary are 1 and N. Now, the field components on the nodes just inside the 
boundary exist in the second page through the (N-l)th page. The second page has 
E x (ij,2,t), E y (ij,2,t) and E z (ij,2,t). For the x compoonents, the third page has [E x (2j,3,t), 
E x (i,N-l,3,t), E x (N-lj,3,t), E x (i,2,3,t)], mth (3<m<N-2) page has [E x (2j,m,t), 
E x (i,N-l,m,t), E X (N-1 j,m,t), E x (i,2,m,t)]. The y,z components have the same global 
coordinates. Finally,the (N-l)th page has E x (ij,N-l,t), E y (ij,N-l,t) and E z (ij,N-l,t). 
However, the limit s on the values of ij,k on the nodes just inside the boundary are 2 and 
N-l. Because the each node on the boundary has three field components, the electric field 


vectors for the grid on the boundary can be formed as below 


Exjboundaryif) ~ \E X _\stj>age E x JXndj)age 

E x jfl-\yh_page E x _Nth_page\ 

(4.75) 

Ey _ boundary (f) = \Ey_lst_page Ey_2nd_page 

>•* Eyjjj-\yh_page EyJJth_page\ 

(4.76) 

E z _boundary(f) ~ \_E z _\st_page E z _2nd_page ** 

... E z _{N-\)thj>age E z j^th_page\ 

(4.77) 


where each page field element forms a row vector with the fields at the nodes on the 
boundary and T indicates transposition (the vector is a column vector, but it has been 
written as a transpose of a row vector to save space). In an analogous maimer an electric 
field vector can be formed for the nodes just inside the grid boundary 


E x JustJnsideif) = \E X _2nd_page - 

— E x _(N-l)th__page\ 

(4.78) 

EyJust_inside(f) = [EyjZnd_page — 

- EyjN-\yh_page\ 

(4.79) 

E z Just_insideif) = \EzJlnd_page . 


(4.80) 

Ejustjnsideif) ~ \E X JusrJnsideif) EyJust_ 

inside (f) E z Just_insidei f) ] 

(4.81) 


where each page field element forms a row vector with the fields at the nodes just inside 
the boundary. The electric field on the grid boundary E^^t) can be expressed as a 
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convolution of the electric fields at the layer just inside the grid boundary E justJnside (t) and a 
matrix of impulse responses H(t). 

Ex_boundary(f) = E x (t) * Ejustjnsideif) (4.82) 

EyJ)oundary(f) Hyif) * Ejustjnsideif) (4.83) 

E z _boundary(f) Ez (0 * Ejustjnsideif) (4.84) 

where the 2(N)(N)^(N-2)4(N-1) by 3[2(N-2)(N-2)+(N-4)4(N-3)] matrix of impulse 



Hy(t) and EL,(t) have the same matrix size as H^t). The electric field(Ex, Ey, Ez) at 
a boundary node (identified by a particular value of the index m) can be written as the 
product of the m-th row of the impulse response matrix H(t) (H^t), Hy(t), H z (t)) and the 
column vector of the electric fields just inside the boundary (indentified by index n) 

E x _boundary(tn, 0 = 2 W= \ h x mytif) * Ejust_insi4e(n, f) = (4.86) 

2 n =i Jq Ejustjnsideif 5 "0 ' h x myi(t ~~ X)cfZ 

Ey boundaiyQn, 0 = S„ = ] hy m.nij) * Ejustjnside(n, t) — (4.87) 

S n= i J Q Ejustjnside (U, l) • hy mJ1 (t ~ X)dz 
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Ez boundaryipi, 0 = 2^=] hz mjiif) * Ejust_inside(n, 0 = (4.88) 

2 n=\ Jq Ejust_inside(ft, "0 • Hz mji(t ~ X)dx 

The discretized forms of the above equations, involving samples of electric fields at t=kAt, 
are the sums of discrete convolutions. Since a discrete convolution is itself a sum, the 
result for the boundary nodes will involve double summation. 

M ■ < 4 - 89 > 

Ey_bou*daiy( m ) ~ d^T=l EpLj ' hymjJ (4.90) 

4,™*(«) (4.91) 

The convolution equations express the fields at the edges as weighted sums (the 

weighting coefficents are the values of the discrete boundary impulse responses hj of the 

time histories of the fields "just inside" the grid. In that respect the concept and 

procedure used in 3-D TGT is the same as those used in 2-D TGT. Also, the impulse 

responses h^t) converge to zero relatively fast which effectively reduces the number of 

significant terms in the convolution sums. The convolution sums can thus be truncated 

such that only a certain number of the most recent values of the electric fields just inside 

the boundary are needed. The discrete boundary impulse responses h mjl (kAt) need to be 

determined (and saved) only once, just like in 1-D, 2-D, for a selected grid velocity 

v d =Al/At. The impulse responses can be stored in a rectangular matrix form. A 

particularly suitable matrix form would have [2(N)(N)+(N-2)4(N-1)] rows and 

3[2(N-2)(N-2)+(N-4)4(N-3)]N h columns( because the field components at the nodes just 

inside the boundary consist of Ex, Ey Ez and the number of each elctric field is 
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2(N-2XN-2)+(N-4)4(N-3)]) , where N h denotes the "truncated" length of the discrete 
boundary impulse responses. This form is suitable because the boundary electric fields can 
then be calculated as a product of such a matrix and a column vector constructed of 
3[2(N-2)(N-2)+(N-4)4(N-3)] subvectors of length N h . The subvectors represent the N h 
most recent values of the electric field at nodes just inside the boundary. The choice of N t 
depends on the desired amplitude "resolution". These subvectors are updated at each time 
step by shifting all of their values down, that is, "futher into the past", by one, discarding 
the last value, and updating the first subvector element with the most recent field value 
calculated via standard grid update equations. The double summations are thus efficiently 
executed as matrix-vector multiplications. 

In the 3-D TGT, the important thing is how to determine the discretized impulse 
responses. The discretized impulse responses hj^ are determined in a manner analogous 
to the 1-D, 2-D cases, using the discrete equivalent of the Dirac deha function which we 
will denote as 8\ In order to determine the discrete boundary impulse response for a 
particular input port n, all input ports except the n-th input port are set to zero and the 
Dirac impulse is applied to the n-th input port. Recording all the outputs from t=0 to 

t=N h At provides us with 2(N)(N)+(N-2)4(N-1) discrete boundary inpulses of duration N h . 

E k xJustJnside (n) = ^ (4.92) 

E\ justjnsidM - 0 for p - \...2(N-2) 2 + (N— 4)4(N- 3) and p*n (4.93) 
E k yJustJnside (p) = 0,E k ZJUStJmide (p) = 0 for p = 1...2(AT— 2) 2 + (AT— 4)4(N— 3) (4.94) 

The h, k ^ and hj^ are determined in the same manner as the h* 1 ^ 

E k yJusUnside (n) = b k (4.95) 
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Eyjustjnsideip) = 0 M P = E..2(N - 2) 2 + (N - 4)4(7V-3) and p*n (4.96) 
E k XJUStJmide (p) = 0 XjustjnsiM = 0 for p= l...2(N- if + (N- 4)4(AT- 3) (4.97) 

El JustJnside (f) = ^ (4-98) 

E k zjustjnside (p) = 0 for p= l...2(N-2) 2 + (N-4)4(N-3) and p^n (4.99) 
jEJj*_ we(P) = 0,^, ^(p) = 0 /or p=l...2(N-2f+(N-4)4(N-3) 

(4.100) 

The each process is repeated 2(N-2)(N-2)+(N-4)4(N-3) times and the results for 
the each case are stored in the 2(N)(N)+(N-2)4(N-1) by 3[2(N-2)(N-2)+(N-4)4(N-3)]N h ( 
[E x E y EJ field vector just inside the boundary) matrix. This matrix is then used to 
implement the transparent grid termination via the matrix multiplication by the 
3[2(N-2)(N-2)+(N-4)4(N-3)]N h column vector of the time-histories of the fields just inside 
the boundary. Depending on the values of N h and N, there will be input nodes that are 
sufficiently far from an output node such that the time it takes for an input impulse to 
propagate to the ouput port exceeds N h . In such a case, the contribution of such an input 
node to the output node is known to be zero in advance. This may be used to reduce the 
n umb er of operations in TGT implementation. The process of obtaining the 3-D discrete 
boundary impulse responses is carried out on a "large" grid whose size should be at least 
(N+N h by N+N h by N+N h ), or equivalently, the distance from the TGT boundary to the 
boundary of the"large" grid should be at least N h Al/2. The termination of this "larger" grid 
(in typical applications N»N h and the "larger" grid would be only incrementaly larger) is 
imma terial, since the reflections off its boundary do not arrive at the output nodes (where 
they would have corrupted the impulse responses) before the time-stepping has been 
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terminated. When calculating the discrete boundary impulse responses (DBIR), the nodes 
interior to the TGT boundary need not be updated since all the nodes on the layer just 
inside the TGT boundary (the input ports) are zero for t>0 (at t=0 only one node on the 
layer just inside the TGT boundary has the value of 1). The DBIR calculations thus 
require updates only for the nodes between the TGT boundary and the large grid boundary 
which can result in significant computational time savings if N»N h . The update 
equations used for the region between the "inner" (TGT) and the outer boundaries are the 
free space equations, since this region is assumed to be free space. (The medium between 
the two boundaries can be any homogenous medium extending to infinity, but the most 
practical case is the free space.) 

C. 3-D TGT RESULTS 

The overall procedure for 3-D TGT testing is very similar to that used for 2-D 
TGT testing. The main difference is in the amount of computer memory needed for 3-D 
TGT and " infinit e" grid implementation. Since 3-D FD-TD calculations require much 
more computer memory because the number of nodes in 3-D is proportional to N 3 as 
opposed to N 2 for 2-D (N is the number of cells on each side of the computational 
domain) this limits the grid size and the duration of impulse responses. We have thus 
selected a domain with only N=ll electric cells on each side (electric cell faces are aligned 
with the domain boundaries). Note that even such a modest number of cells per each side 
results in N 3 =1331 electric cells for the domain volume. Since each cell has three field 
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components, the total number of electric field components is 3N 3 =3993. The number of 
magnetic field components is about 3(N-1) 3 =3000. The total number of unknowns is thus 
about 7000. 3-D TGT is implemented by updating the electric field components on the 
TGT boundary based on the time histories of electric field components at nodes just inside 
the TGT boundary. The source is a delta impulse of electric field E field at the center of 
the grid, with unit-amplitude components E x , E y and E z . The source creates an outgoing 
spherical wave. If the grid and the TGT were perfect the power within the TGT will be 
constant until the outgoing wave has reached the nodes in the centers of the six TGT 
boundary planes (because these nodes are the closest to the source located at the grid 
center. From then on the power within the TGT should decrease until the outgoing wave 
has reached the comers of the domain (the nodes furthest away from the source at the 
domain center). After this time the power within the grid, for a perfect grid and TGT, 
should be zero. However, since the grid is nor perfect, there will be some residual power 
within the TGT boundary even for an infinite grid. This is the result of the impulse 
"spreading" as it propagates through the FD-TD grid. Furthermore, the residual power 
will include the power "reflected" off an imperfect TGT. Again, just like in 2-D, the 
residual power can be calculated as a function of time for the two cases infinite grid and 
TGT. The relative difference of the two expressed in dB can be thought of as a measure 
of performance of the TGT. This quantity can be plotted vs. time (with the time origin 
defined as the tim e when the outgoing wave first reaches the TGT) for various impulse 
response durations. Unfortunately the results for 3-D TGT were not satisfactory. The 
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TGT tests showed "late-time instability", that is the power within the grid would at first 
decrease (as it should) but later will start to increase without bounds, which is incorrect. 
A sample result is shown below (what is the duration of the impulse response for this?) 



Figure 4.14. Residual Power Within TGT Showing Late-Iime Instability. 

The non-physical oscillatory behaviour of the solution is also evident from the plot 
of the electric field at a node inside the TGT. The field at first converges to zero as it 
should, but after about 200 steps it starts to diverge in an oscillatory manner. 
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Figure 4.15. Field Strength at a Node observation point within boundary. 


For an infinite boundary, the residual power within the TGT boundary decreases to 
zero as it should. The late-time instability could not be observed because the available 
PC's did not have enough memory to accomodate the size of the required "infinite" grid. 
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Figure 4.16. Residual Power within TGT for an Infinite Boundary. 
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V. SUMMARY AND CONCLUSION 


A. SUMMARY 

In this thesis we have tested the performance of the TGT in 1-D, 2-D and 3-D, 
based on the residual reflection off the TGT boundary. The power within the TGT 
boundary was calculated at each time step for an "infinite" grid (that is a grid so large that 
the reflection off its boundary do not reach the nodes within the TGT boundary before the 
calculations are completed) and for a much smaller grid with the TGT. For 1-D there were 
no reflections off the TGT boundary, that is, the outgoing wave was "absorbed" perfectly 
by the TGT. Note that in 1-D there is no impulse spreading. Therefore, in 1-D, the grid 
and the TGT can be "perfect". This is not the case in 2-D however. Neither the grid nor 
the TGT can be perfect, that is there is some pulse spreading and some "reflection" off the 
TGT. However, the power reflected off the TGT is less than 1% of the residual power 
due to pulse spreading. Therefore, although the 2-D TGT is not perfect, its performance 
is excellent since the reflections off the TGT are well below the level of the "distortion" 
introduced by the grid itself Furthermore, the power reflected off the TGT boundary 
depends on the duration of the truncated impulse responses: truncating impulse responses 
at a later time (increasing the impulse response duration) reduces the reflections off the 
TGT. Finally, the 3-D TGT implementation was not successful, because of late-time 
instability problems. 

B. CONCLUSIONS AND RECOMMENDATIONS 

The results of this study show that TGT concept works for 1-D and 2-D FD-TD 
electromagnetic field problems. The 3-D TGT implementation was the most difficult, 
because of the computer memory requirement. Unfortunately, it was not successful due to 
late-time instabilities. More research is needed to determine and remove the causes of 
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late-time instabilities in 3-D. A previously completed thesis [Wells] has demonstrated the 
effectiveness of TGT for 2-D static problems, while this thesis has demonstrated the TGT 
effectiveness for 2-D time-domain problems. Another possibility for further research is to 
implement the TGT in the frequency domain, where the time domain convolutions would 
be replaced by sums of products of complex numbers. 
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APPENDIX 


A. DBIR.M 


%%%%%%%%%%%%%%%%%%%%%%%%*%****%%%***%%********************* 


%%%%% Program : DBIR.M ***** 
%%%%% This Program Generates TGT Coefficients for 2-D %%%%% 
%%%%% TMz fields (Ez, Hx, Hy) ***** 
%%%%% NOTE: All H-fields have.been multiplied by ZO I %%%%% 


%%%%%%%%%%%%%%%%%%%%%%%*%%****%%*%*******%***************** 


clear all 
clear, flops(0); 
tic 


N = 101; 

Nin = 21; 

Nout = Nin + 2; 


h_dur - 40; 
vr = l/sqrt(2); 
rowsh = 0.5*(N - Nin); 
colsh = 0.5*(N - Nin); 
irt = rowsh + 1; 
irb » irt + Nin - 1; 
icl = colsh + 1; 
icr = icl + Nin - 1; 
nodes__out = 4*(Nin+l); 
nodes in = 4*(Nin-1); 


% duration of impulse responses 
% grid "speed" 

% Row "shift" for the "inside" matrix 
% Column "shift" for the "inside" matrix 
% Top row for the "inside" matrix 
% Bottom row for the "inside" matrix 
% Leftmost column for the "inside" matrix 
% Rightmost column for the "inside" matrix 
% Number of nodes on the boundary 
% Number of nodes just inside the boundary 


%%%%% INITIALIZE ARRAYS %%**** 

Ez = zeros(N,N); Hx = zeros (N-2 ,N-1) ; Hy = zeros(N-l,N-2) ; 

Edbir = zeros (nodes_out, h_dur) ; Dbir = zeros (nodes_out, nodes_in*h__dur) 
row pointer = zeros (1 ,nodes_in) ; col_jpointer = zeros (1, nodes^in) ; 
col_start = zeros (1, nodes_in) ; col_stop - zeros {1, nodes__in) ; 

for m = l:l:nodes_in % LOOP OVER INNER GRID NODES 

if (m <= Nin) 

row__pointer(m) = irt; col_pointer(m) = icl + m - 1; * 

elseif (m > Nin) & (m <= 2*Nin-l) 

row_pointer (m) = irt + m - Nin; col_pointer (m) * icr; * 

elseif (m > 2*Nin-l) & (m <= 3*Nin-2) 

row_pointer (m) = irb; col_pointer (m) = icr - (m-(2*Nin-l)) ; % 

elseif (m > 3*Nin-2) 

row_pointer (m) = irb - (m-(3*Nin-2) ) ; col_pointer (m) = icl; % 

end 

col_start(m) = (m-1)*h_dur+l; col_stop(m) = col_start(m) + h_dur - 

end 

pack 


%**%*************************************** 
%%%%%%%%%% DISCRETE BOUNDARY IMPULSE RESPONSE LOOP l %%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%*%%**************% 


% LOOP OVER INNER GRID NODES as sources 
for m = l:l:nodes_in 

Ez * zeros (N,N) ; Hx = zeros (N-2 ,N-1) ; Hy = zeros(N-l,N-2}; 


TOP SIDE 
RIGHT SIDE 
BOTTOM SIDE 
LEFT SIDE 
1 ; 
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Ez (row_pointer (m) , col_jpointer (m) ) - l; 


% Delta pulse excitation along the. top 
half of the boundary 

%%%%%% TIME LOOP %%%%%% 
for k = 1:1:h_dur 

Etop = (Ez(irt-1,icl-1:icr+1))'; 

Eright = Ez(irt:irb+l,icr+1); 

Ebot = (fliplr(Ez(irb+1,icl-1:icr)))'; 

Eleft = flipud(Ez(irt:irb,icl-l)); 

Edbir(: , k) = [Etop; Eright; Ebot; Eleft]; 

Hx = Hx + vr*( Ez(2:N-1,1:N-1) - Ez(2:N-1,2:N)); 

Hy = Hy + vr*(-Ez(1:N-1,2:N-1) + Ez(2:N,2:N~1)); 

Ez(2:N-1,2:N-1)=Ez(2 :N-1,2 :N-1)+vr*(-Hy(1:N-2, :)+Hy(2:N-1, :) 

+Hx{:,l:N-2)-Hx(:,2:N-1)); 

Ez(irt:irb,icl:icr) = zeros(Nin,Nin); % Reset Ez inside 
end 

%%%%%% EKL OF TIME LOOP %%%%%% 

Dbir(:, col_start (m) : col__stop (m) ) - Edbir; 

end 

%%%%%% END OF SOURCE LOOP %%%%%% 

MFlops = flops/1000000 

toe 

save dbir__40m Dbir 

save par_40m Nin h_dur vr 


B. TGT TM.M 




%%%%%% Program : TGTJTM.M %%%%%% 
%%%%%% FD-TD for 2D, TMz fields (Ez, Hx, Hy), cartesian cordinates %%%%%% 
%%%%%% NOTE: H-fields have been multiplied by ZO ! %%%%%% 
%%%%%% Dirichlet b.c. or OPEN boundary simulated via DBIR %%%%%% 




load dbir_4 Om 
load par_40m 

%%%%% INPUT PARAMETERS %%%%% 

N = Nin + 2; 

nodes out = 4*(Nin+l); % Number of nodes on the boundary 

nodes_in = 4*(Nin-l); % Number of nodes just inside the boundary 

% Time Samples 
N_samples = 200; 

%%%%% ARRAYS INITIALIZED %%%%% 

Ez = zeros(N); Hx * zeros(N-2,N-1); Hy = zeros(N-1,N-2); 

Ebound = zeros (nodes_in*h_dur, 1) ; E_dbir = zeros (nodes_out, 1) ; 
time_pointer - zeros(nodes__in,1); 

%%%%% BOUNDARY POINTERS FOR THE TIME HISTORIES OF NODES JTJST INSIDE %%%%% 

for i = l:nodes_in, time^pointer(i) = (i-l)*h_dur + 1; end 

%%%%% BOUNDARY POINTERS FOR SPIRAL TO RASTER LABELING CONVERSION %%*%% 
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%%%%% (can be done only once and stored!!!) %%%%% 

for n = 1:1:nodes_out 
if n <= N 

row_pointer (n) = 1; coljpointer (n) = n; 

elseif (n > N) & (n <= 2*N-1) 

row_pointer (n) = n - N + 1; coljpointer(n) = N; 
elseif (n > 2*N-1) & (n <= 3*N-2) 

row_j?ointer (n) = N; col_pointer (n) = 3*N - 1 - n; 

elseif (n > 3*N-2) 

row_pointer(n) «= 4*N - 2 - n; col_pointer (n) = 1; 
end 

end 



%%%% Computation of Residual Power for grid %%%% 

%%%%• since wavefront touching the tgt boundary %%%% 

if (k>=ceil(N/2) & k<=ceil (N/2)+140) 

m=m+l; 


133 



s (m) = 0 ; 
for a=l:N 
for b=l:N 

s(m)-s(m)+Ez(a,b)^2; 
end 

end 


end 


end 

save s_tgt s 
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